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AEROFOILS OF MAXIMUM THICKNESS RATIO FOR 
A GIVEN MAXIMUM PRESSURE COEFFICIENT 


By A. R. MANWELL (University College, Southampton) 
[Received 3 October 1947. Revised 2 December 1947] 


SUMMARY 

A relation is found which determines the symmetrical aerofoil which in streaming 
incompressible or subsonic compressible flow has the greatest thickness ratio for 
& given maximum pressure coefficient. The solution in the incompressible case 
coincides with that of a free streamline problem due to Riabouchinsky. The 
Kaérman-Tsien solution is derived in a very similar form. It is also shown how in 
principle one experiment with an electric potential apparatus would solve the 
problem for any pressure-density law. 

The method of deriving the relation can be extended to various similar problems, 
including axially symmetric flow, and the solution possesses the further property 
of maximizing the area and other functions of the contour of the aerofoil. 


I, Introduction 
Tue problem to be solved is that of finding the shape of the symmetrical 
aerofoil which in unlimited streaming flow has the following properties: 

For a given maximum pressure (or velocity) peak on the boundary 

(a) the maximum thickness to chord ratio, 
(b) the area to (chord)? ratio 
are as large as possible. 

Lockwood Taylor (ref. 1) argued that the most efficient method of 
obtaining lift would be to make the suction, and so the velocity, constant 
over the boundary. Using thin aerofoil theory he calculated the approxi- 
mate form of both cambered and symmetrical aerofoils having constant 
velocity portions. In an unpublished paper, L. G. Whitehead in 1942+ 
has worked out a number of similar aerofoils of finite thickness having 
constant velocity arcs. He observes in his paper that a free streamline 
problem solved by Riabouchinsky (ref. 2) has, of all other examples known 
to him, the best thickness ratio in the sense defined earlier in this para- 
graph. The present paper is, so far as the author is aware, the first to 
attempt to formulate and discuss the problem precisely. Apart from 
giving a proof for incompressible flow, the method used can be readily 


generalized to subsonic flow, by using a well-known theorem of analysis. 


tI am indebted to the Secretary of the Aeronautical Research Committee for the 
pportunity of reading Whitehead’s paper, entitled Minimum Velocity Aerofoil, 6121, 
Ae 2073 (dated 15 August 1942). This hitherto confidential paper will shortly appear as 
#Report and Memorandum of the Aeronautical Research Committee. 

5092.4 Bb 
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The conjecture made by Whitehead (loc. cit.) that Riabouchinsky’s solu- 
tion is very near to the true solution for compressible flow is now replaced 
by the proof that a certain shape defined in the hodograph plane in a 
very similar fashion is truly the best possible and it reduces to the solution 
of ref. 2 for incompressible flow. Some numerical comparisons are made 
for incompressible flow, but, unknown to the present author, these had 
already been carried out in greater detail by Whitehead. They are given 
for the sake of completeness (Fig. 2). 


II. Two-dimensional incompressible flow 
Since the aerofoils considered are symmetrical, only the upper half- 
plane need be considered and the flow will be taken along the x direction 
with stagnation points at each end on the real axis. To make the problem 
definite the aerofoil must be restricted to lie between two lines, say x = +1. 
Let an infinitesimal bulge be made at a point D on the boundary of the 
aerofoil. In physical terms the effect of a bulge is like a doublet with 
its axis tangential to the boundary at D, and with the old boundary and 
the real axis as given streamlines. It is evident that the flow lines due 
to this disturbance are such as to decrease the velocity on the original 
boundaries. 
More precisely, mapping the upper half of the physical plane (of 
z = x+iy) on the upper half of the f-plane so that the aerofoil goes into 
the strip BC, Fig. 1, and taking the image of D in the ¢-plane as origin, let 
a new plane be defined by: 
& (1 —e?/#?)(1+-0%e?/t?), (1.1) 
where « is infinitesimal. 
The semicircular bulge t = ee? maps back into the ¢-plane and, in 
general, into the z-plane as a bulge without stagnation points, and which 
can be made as flat or as sharp as is desired by choosing 0 < a? < 1. 
Let the complex potential be 
w= UL (1.2) 
for the undisturbed flow and 
w, = U(t+e/t) (1.3) 
for the disturbed flow. Then, if Q and Q, denote the corresponding com- 
plex velocities 


Q _ dw, _ dw dw, _ o_ 1—e*/? _ 
1 dz dz dw = * (1—e?/#®)(1 +. a%e?/t?)’ 
na 1.4 
and hence 12,| = |Q| ral: (1.4) 
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Except at the bulge, ¢ is real on the boundary and therefore 
@| < |@|; (1.5) 
at the bulge we have 
l 
V/{(1+a?)?— 2a? cos 28} 


Q; @ (1.6) 


9 


where the factor varies between 1/(1-+-a?) at the ends of the bulge and 
\/(1—o?) at the middle of the bulge. This infinitesimal variation of the 
boundary has the effect of 

(a) increasing the area of the aerofoil; 

(b) increasing, or leaving unchanged, the maximum thickness; 

(c) decreasing the velocity except at the bulge. 

If the aerofoil has either of the required properties, it follows that at 
D either 

(i) the velocity was originally the maximum, say M, so that at the 

bulge |Q,| > M however small a?, or 

(ii) the chord was the maximum and D lies on x = +. 
It is concluded, therefore, that an aerofoil possessing either property (a) 
or (b) stated in § I must be such that at every point the velocity or chord 
isthe maximum. This property defines an aerofoil in the hodograph plane. 


III. Two-dimensional compressible flow 

The same result will follow in compressible flow if it can be shown that 
the effect of a bulge is that of a doublet. That this is so in the subsonic 
case is seen as follows: 


The equations of two-dimensional flow are: 


a. (2.1) 
Cy ea 
< (pu) +—(pv) = 0 (2.2) 
OX Yy 
iO, sins ue ‘ 
+ 43(u*+v*) = constant, (2.3) 
J P 
where p = p(p) (2.4) 
and wu, v are velocity components in the z- and y-directions respectively. 
Let pu —y, and pv=y, (2.5) 


where the suffixes denote differentiation with respect to x and y; (2.2) is 
then satisfied and (2.1) gives 


(2.6) 
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Expanding and using (2.3) and (2.4), 


1 dp 
q dq 
e+ yy — —2 [hee Wey ey Hin] = 9% (2.7 
(99) 7, (00) 
where g? = u?+-2?, 
Defining the speed of sound ¢ by = = c*, then equation (2.7) becomes, 
ap 4 ry 
using (2.3), (2.4), and (2.5), HY J Ay | 


(te2+Pyy)[p? (q 2 ¢2) )]—W2 b+ I 2h. by Pry the by) = 0. (2.8) 


By (2.3), (2.4), (2.5), p and q are functions of ¥2+-42 only, and the equation 
(2.8) is of elliptic type if 


[e*(q°—c?) pz] p°(q?—c*) pa] oy > 0, (2.9) 
or pee Re ] >9, 
i.e. if q* < c*, (2.10) 


If the bulge is now made on the boundary as in the incompressible case 

and « is infinitesimal the new stream function may be written 
wy = +n, (2.11) 
where 6 is a small constant. 

Except for the stagnation points in the original flow the first derivatives 
of 5y are always small compared with those of %. Then in a region limited 
by the aerofoil with bulge, small circles enclosing the stagnation points, 
and the real axis, 7 satisfies an equation of the form 


A nye t+ 2Bygy+ OC nyt Engt+ Fr, = 9. (2.12) 


Here A, B, C, EH and F are known analytic functions of (2, y) because of the 
known analytic solution of (2.8) representing the original subsonic flow 
and since AC > B?, by (2.8), (2.9). Also 7 is zero on the real axis and 
the boundary of the aerofoil and, since the original streamlines cross the 
bulge and the circles enclosing the stagnation points, is of one sign, say 
positive, on these parts of the boundary. Finally, on a sufficiently large 
semicircle 7 tends to zero independently of the other small quantities in 
the problem. By a well-known theorem (ref. 3) the solution of (2.12) has 
no maximum or minimum in a region where it is regular, and it follows 
that 7 is positive in the region. Therefore on those parts of the boundary 
where 7 vanishes the normal derivative (@7/én) > 0. By considering the 
behaviour near the ends of the bulge it is seen that the new velocity 
(@b,/On) < (@/@n) at all points on the boundary. 
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IV. Various extensions of the results 


The first extension to be made is to axially symmetric subsonic flow. 
This may be visualized as a pseudo-two-dimensional flow in the (2, y)- 
plane with y > 0. The existence of the disturbing doublet flow has been 
seen to depend only on the fact that the differential equation for y is 
of elliptic type. 

The differential equation for the stream function is 


=(- *\+aloa) = 0. (3.1) 
éx\py dx} © dy\py oy 


The second-degree terms reduce to 
(us. Pyy)p?q?(1 -c?/q?)—(1 Y?) (UE eet he by Pry t+ by): (3.2) 
The equation is of elliptic type if 


9 


[ p?(q? c*) of” y* || p?(q? -c*?)— y?| > 


_ vey « 
” 
which reduces to gf <c. (3.3) 


The boundary conditions for flow past a closed body of revolution are 
the same as for symmetrical two-dimensional flow and so the results of 
§$II and III are again true, namely that at every point the velocity is 
the maximum or the point lies on x +], 

The next extension which can be made to all the preceding cases is that 
involving the presence of symmetrical disturbances in the stream, such 
as vortices or fixed bodies. For example, the relation of maximum chord 
or velocity applies to an aerofoil body symmetrically placed in a tunnel. 

Finally, the method of proof shows that any solution obtained maxi- 
mizes not only the thickness and area ratios but all functions which are 
additive with area, such as the moments of any order about the z-axis. 


V. Solution for incompressible flow 


The incompressible case is essentially the solution given by D. Ria- 
bouchinsky (ref. 2) for a flow with free streamlines. 
Let Q dw/dz; then the complex potential corresponding to equa- 
tions (3) and (6) of that paper is given by 
0+ Q-)/(u—u 
( (+e uk oll Rais (4.1) 
(1+[(Q+Q)(u-*—u) P} 
where iw is the complex velocity at infinity. The hodograph plane with 
streamlines is sketched in Fig. 1 and it is readily verified that equation 
(4.1) takes real values on the boundary. 
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Corresponding to equations (7) of ref. 2 the boundary is given in 
parametric form as 
a sin 6 cos 0 cos?6 dé 
~— (e2-+-c0s?6)# (e2-+ cos?6)?’ 
0 


1 

1—q?) dc , , ce 

iO) = Set | ream gtamT Ltey t+ (et cost) 
0 

(4.2) 














¥ B'B c 


Fic. 1. (i) The physical plane, (ii) the hodograph plane. 
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Constant velocity, AA ; ellipses, BB; circular ares, CC ; Rankine ovals, DD. 
where e = }(u-1—u). The first term of y(@) corresponds to a straight 
portion at the nose, and @ is the coordinate of a point on the unit circle 
in the hodograph plane. 
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AEROFOILS OF MAXIMUM THICKNESS RATIO 371 
Finally, the thickness ratio, cf. equation (14) of ref. 2, is given by Y/X, 
where 
X = {1/(1+e*)!}{(1+-e*) E[1/(1+e?)#]—e?K[1/(1+-e?)*]}, 


1 


2 [ (1—q?) dq : . 
- ms —e2/(1 + e2)8 43 
| [(1+-u*q?)(1+-¢, wp’ e/(-+e’) (4.3) 





and LE, K are standard complete elliptic integrals. 
As e > 0 the peak velocity approximates to u(1+e) and X ~ 1, Y we. 
It follows that for all symmetrical aerofoils we have the following: 
For a velocity peak of u(1+-e) where e is infinitesimal and U is the stream 
velocity the best thickness ratio possible is e. 


VI. Calculation of the profile for subsonic flow: in particular the 
Karman-Tsien case 
With the notation of § IT, let 
vt ™ v= dy, (5.1) 
and let g and @ be the magnitude and inclination of the velocity to the 
y-direction. It is well known (refs. 4, 5) that ¢ and %& satisfy a pair of 
linear equations with g, 6 as independent variables. 
Thus it can be easily verified from (2.1) to (2.4) that 
| od d l Ous ~ 
2 ean, (5.2) 
q oq dq \pq} 08 
/ Ou ] od rk Q 
pyoqg  g? oO a 
The consideration of many problems is aided by the observation that 
(5.2) and (5.3) can be identified with the equations of incompressible flow 
ina thin stratum in the plane of (¢,@). Thus, (5.2), (5.3) are identical with 


Cd a L (5.4) 
or hr 00 
l od l + (5.5) 
r 06 h or 
where h and 7 depend on gq only. 
Comparison of (5.2), (5.3) with (5.4), (5.5) gives 
h = p(l—q?/e?)-+, (5.6) 
l 
mee [“0-¢ ci (5.7) 
* Pp 


The usual pressure-density law can be written 


p = kp’+constant, (5.8) 
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from which it follows that 











— ‘ 1(y—1) —} 
h/h, = ( —~—¢ <é) (1 ¢’ g—*0') (5.9) 
t—1\4/t-+ vp\"2v8 . 
Tlf. = | —— o.1 
. (=) fae aie 
1—4(y—1)q?/c?2 —l 
where i? — <a. Dace § = i, (5.11) 
1—}(y+1)q?/cG y+l 
and where the suffix 0 refers to the value for g = 0. 
For an isothermal flow (y = 1), we have 
pf —,/(1—g2/c2)} 
rity = gio Aaa 8 os (5.12) 
1+./(1—q*/c*) 
h/hy = (1—q?/c?)-te-4°e*, (5.13) 
whilst the Karman-Tsien approximation y = —1 gives (cf. ref. 6, 7) 
h/h, = 1, (5.14) 
r/ro = (q/c)/{1+.(1+¢?/c?)}. (5.15) 


For a general adiabatic law the flow in the plane of (r,@), which might 
be called the modified hodograph plane, is seen to be identical with the 
incompressible flow in a stratum of depth A(r), where h(r) is defined by 
(5.9), (5.10), (5.11). 

The representation is true only for the subsonic case, which is just the 
case for which in § III it was shown that the required aerofoils had the 
simple defining property of maximum velocity or maximum chord at 
each point. Then, 

dd = ¢,dx+¢4,dy = —qsin@dx—q cos 6 dy, 
dis = b,dx+i,dy = pq cos 0@dx—pqsin 6 dy. (5.16) 
Along a streamline therefore ds = 0 and 
dx = —dd¢dsin6/q, 
dy = —dd¢cos6/q. (5.17) 


If the potential satisfying equations (5.4) and (5.5) is known, these equa- 
tions determine the physical boundary. Since ¢ and % satisfy the same 
boundary conditions in the subsonic case as for incompressible flow, the 
problem is in principle solved, e.g. by constructing the flow in an electrolyte. 
In the Kérman-—Tsien case an explicit solution is easily derived. That 
this approximation is valid follows from the fact that h(q) varies little over 
a wide range in both adiabatic and isothermal flow, though, in both 
cases, it eventually becomes infinite (see Table). 
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Table of the functions h(q) and r(q) 























{diabatic flow Isothermal flow 
Co 9/9 h r q/e | h r 
o I*0o ° ° | I'oo-) | ° 
"22 I*00 *20 “2 | 100 *20 
4 14 ‘OI *40 4 I°or +38 
66 1:06 *56 6 1°03 “54 
8 88 1°38 *69 8 1°20 *67 
9 50° ‘TI 9 1°56 | — 
I I x 71 1°00 io 8) | 73 
q al fluid velox ity ‘ 
fo peed of sound at stagnation point. 
1 itical speed = cy | ——. 
i 
Consider the function 
W—o+i¥ eee Fh (5.18) 
> > 4 7\|2\? 5 
Ja+[(P+P>)/(V>— VF} 
where V r(u)/r(M), Fr re9/r(M) = pet, (5.19) 


*) 


and M is less than co,/| : ) the critical velocity, and r(q) is defined by 


equation (5.15). Then W satisfies the boundary conditions in the modi- 
fied hodograph plane, and so the exact solution on the Karman-Tsien 
approximation has been found. 

, an CO . 2F io 10e¥ m 

Now Le = —— —4+— _., (5.20) 

dP -° oe P 00° P 00 

By integrating along the contour corresponding to the boundary of the 
aerofoil the following is found, corresponding to § V, equations (4.2). Write 
E = 4(V-1—V); then the straight portion BC (Fig. 1) is of length 


710 td (1 p*)p dp : : 
AH® 72 9/179\ 18 7 ) say. (5.21 
| (1 ps *)\(1 + p* J 2) |8q(p) ne ) 


and putting in 





= ; i sin*6—sin 6 cos 6 
X(0)—Y(0) = E? E —_____—_____ 5.22 
ii ical (£?+-cos?@)# ; \ ) 
the maximum thickness-to-chord ratio is Y,/X, where 
Y, = Y(4m7)+, (5.23) 
X, = X(}z). (5.24) 


For infinitesimal thickness ratio 


X,2>1 and Y¥,>E#. 
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Let the peak velocity M = U(1+e) and therefore 
, r(u) 1 dr 
v= ~ 1—Me|- — 
r(M) (: a 
7), Vw 1—e(1—M2/C%,)t ~ 1—E, (5.25) 
where 0%, = (dp/dp), for the Karman-—Tsien relation between p and p. 





and, using (5. 


~ 


For Karman-Tsien flow the following holds: For a velocity peak of 
U(1+e), where U is the velocity of the stream and e is small, the best thickness 
ratio possible is E = e(1—M?/C%,)}. 


Conclusion 

Figure 2 shows the curve of velocity peak against maximum thickness 
ratio for a number of well-known profiles in the case of incompressible 
flow. The numerical advantage over the ellipse, for example, is slight at 
normal thicknesses. However, since in general the effect of increasing the 
Mach number is to increase the ‘effective’ thickness ratio, it is quite likely 
that the advantage is really rather greater. A complete numerical in- 
vestigation of the solutions for the adiabatic law with higher Mach 
numbers would show the highest possible Mach number attainable} with 
(i) given thickness/chord, (ii) given area/(chord)? ratios. The calculation 
of &% is simpler than for the ellipse or other aerofoil shapes since the hodo- 
graph boundary is known in advance. This property in fact suggested 
the solution before it was realized that it had other special characteristics. 
As an approximation it is suggested that y% be represented by a sum of 
terms, each being the imaginary part of a function of form (5.18) where 
r(q) is any increasing function in the range 0 <q < M; the functions 
being chosen to make the mean square of the error, after substitution in 
the differential equation for ys, as small as possible. 

Finally, I wish to acknowledge the help and encouragement of Dr. 
W.H. J. Fuchs in many discussions of the problem of this paper. He has 
pointed out that a general proof, in the hodograph plane, of the properties (a), 
(b) of § I is difficult because of the possibility of non-simple boundaries. 
A proof of (a) using the hodograph plane, which was found before that 
of §§ II, III, is therefore omitted as being less general than that given. 
The work itself was carried out whilst I was a member of the staff of the 
University College, Swansea. 


+ That is to say, with purely subsonic flow. 
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ON THE TOTAL REFLECTION OF PLANE WAVES 
By F. G. FRIEDLANDER 
(Department of Mathematics, The University, Manchester) 


[Received 17 November 1947] 


SUMMARY 


The reflection and refraction of transverse plane waves at an interface paralle| 
to the direction of polarization is considered when the incident wave is of arbitrary 
shape and the angle of incidence exceeds the critical angle. It is shown that the 
solution of this problem depends on the determination of a plane harmonic function 
h(€, n) satisfying the condition 


(3),_.-A5)__, = 26), 


where A is a known constant and f’(é) a given function. By using the half-plane 
analogue of Poisson’s formula, A(é, 7) can be expressed in terms of f’(E). 

The results show that the reflected and transmitted disturbances exist every- 
where at all times even when the incident wave has a well-defined front, and that 
the transmitted disturbance due to an incident simple pulse is of the order of the 
reciprocal of the distance from the interface, when this distance is large. 

It is pointed out that the same analysis can be applied to the treatment of the 
total reflection of electromagnetic waves of arbitrary shape. 

Finally, the propagation of waves of arbitrary shape over the surface of a semi- 
infinite elastic solid is considered and shown to be possible when the velocity of 
propagation is that of Rayleigh waves. 


THE usual treatment of the oblique reflection and refraction of plane waves 
which satisfy the wave equation can be worked out in terms of infinite 
sinusoidal wave trains or in terms of plane waves of arbitrary form. But 
when the angle of incidence exceeds the critical angle, the analysis, when 
applied to arbitrary wave forms, breaks down because the assumption 
that the reflected and refracted waves are proportional to the incident 
wave is then no longer valid. The object of this note is to set out the 
solution of the reflection problem in this case. 

To begin with a definite example, we shall consider the refraction of 
transverse plane waves in an elastic solid at an interface parallel to the 
direction of polarization. Let the propagation vector of the incident wave 
lie in the zy-plane. Then the displacement, w, is at right angles to the 
ay-plane, and depends only on 2, y, and the time ¢t. Taking the z-axis 
to be the refracting interface, let a’, u’ denote the velocity of propaga- 
tion and the modulus of rigidity respectively for y < 0, and a, pu be the 
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corresponding quantities for y > 0. The displacement satisfies the wave 


tions os a 
equa Cw . ew 1 dw 
—+— =a «Cy <9), (1a) 
ox* —Oy* a’? at? 








ew ew 1 ow 


dx? ' ay? =a At? 








(y > 0). (1b) 


The only non-vanishing stress components are given by 


ow ow 


T... = p—; T,.=p— (y>09); 
™ Ox - oy ty 
sip . i= Ph an y <0). 
zz B ox yz ey (y ) 


At the interface, y = 0, the displacement itself and the stress component 
T,, must be continuous, so that 


. Ow . ,Ow : : 
lim ws ) lim |p =). lim w= lim w. (2) 
y—>+0 oy y—>—0 oy y>+0 y—>—0 
We can take the incident plane wave in the form 
xsiny’+y cosy’ 
wy = ft penn (y <0), (3) 
a 


where y’ is the angle of incidence. Then we can assume 


w=wtw, (y <9); w=w, (y>O0). (4) 





The reflected wave w, satisfies (1a) and the transmitted wave w, satisfies 
(1b). The conditions (2) become 


, [ew ow Ow * 
Wot, = Wz; pw (—®+—t)=p—* (y= 0). (5) 
wy yy cy 
At y = 0, both wy and @w,/éy depend only on 
xsiny’ 
é=t—-—". (6) 
a 


This suggests that a solution can be found by assuming w, and w, to 
depend only on y and é. For the reflected wave w, this leads immediately 


to the form silt cated 
cs — -08 “ 
wy a(t— SUEY), (7) 


a 
where g is an as yet undetermined function. It is also a plane wave, and 
its propagation vector is obtained from that of the incident wave by 
reflection at y = 0. Now the velocity with which the disturbance due to 
the incident wave travels along the interface is 


, 


C= = 7° (8) 
sin y 
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If c >a, a solution of (1b) depending on y and é only is simply a plane 
wave and the usual elementary theory of refraction results. But when 
the angle of incidence exceeds the critical angle, sin-'(a’/a) (which can 
only happen when a’ < a), the elementary treatment breaks down. This 
is not surprising; for in this case c < a, so that the transmitted movement 
caused by the incident wave travels ahead of the incident wave along the 
interface. As the incident disturbance comes on the interface from in- 
finity, it follows that the transmitted wave must exist for all time and 
give rise to a reflected wave, also existing for all time, even when the 
incident wave has a well-defined front. The elementary theory is incapable 
of representing such conditions. 

Considering now the case ¢c < a, we can take, in accordance with our 


Ws = Wt—2 , vI(e -=}}= = h(€, 7), say. (9) 


Substituting in (16) it follows that 
eh 2 *h _ 
ag ont am 
Hence f is a plane harmonic function of £, y, defined for 7 > 0. The 


condition which h must satisfy at 1 = 0 is obtained by substituting (9) 
and (7) into (5), which gives 


f(E)+9(€) = h(E, 0) (10) 


COBY’ gpg ney 1_ 1) (eh 
od f O'eV10) = Ja) 


Differentiating (10) and eliminating g’(£), it follows that 


oh oh aie 
(2), o Mea) o = PO 


a’ 1 1 /(a? sin?y’ —a? 
where A= — = ra sian _ by ( Y ) 
pe cosy ec @ peacosy 
As the angle of incidence, y’, increases from the critical angle to }7, the 
parameter A increases from 0 to oo. 
Assume that the plane harmonic function A has the form 


assumption, 











(11) 








h(&, 7) Lam f He $(a)+% ce Ao) ay (C = €+in, » > 0). (12) 


Differentiating, it follows (using the Cauchy-Riemann equations) that 
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This can be transformed by partial integration. The integrated terms can 


be assumed to vanish, so that 


oh oh | G(x) +h (a) gy, 
0€ C7 «|, 


Separating real and imaginary terms we find therefore that, for » > 0, 


ch . (ee EW) gy, (12a) 
c& ot (a—£)?+-7° 
ch L fi E)h ae 
a 5 )p . 126 


Now if x(a) is a function continuous at « = € and satisfying suitable 
conditions in the neighbourhood of €, it is known (1) that 


x 





; ' (cx) da 
lim 7 x ). 
co (13) 
jm 1 f (Alaa) de _ 1 pf x00) gg, 
nom J (a—£)?+7? w J a—g 


where the prefix P denotes that the last integral must be interpreted as 
a principal value, and must exist in this sense. Applying these results to 
12a,6) it follows that 


(=) $'()+-P je de, 
7=0 -€ 








0€ Tr J a- 
pare (14) 
(=) wie t=P | 2% da. 
On] n=0 7 x—& 
Substituting in (11), 
dh’ (€) +A’ (E) 4 Lp fe ( a ie - 2f'(&). 
—£ 
Hence the function A given by (12) will satisfy (11) if we take 
YF: S 2A F 
p(s) ia I OS): b(€) pi (é). 
| 2 f 9+\e—6) : 
whence h(é, 7) , | > the -* f(x) da. (15) 
oy +r*) J (a— Eft! 
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To complete the formal solution of the reflection problem we need the 
value of h for 7 = 0. Using (13) again, it is found that 





limbén) = 2 f4+—2— Pf £2 ae 
7-0 1+ A? a(1+A?) a—é 


—c 


and hence, substituting in (10), the reflected wave is obtained as 


1—2 2d * f(x) 
Trwl O+ — 


—-o 


g(é) = 








(16) 


The analysis by which we have obtained (15) is purely formal, but it 
is not difficult to state conditions on f(é) which justify the steps by which 
(15) is derived. For example, we can assume f(€) to vanish for large values 
of |€| (which involves little loss of generality in a physical problem) and 
to possess a continuous derivative satisfying a Lipschitz condition: it then 
follows that h, as given by (15), is a plane harmonic function and satisfies 
(11). 

Alternatively, we may simply substitute a given f(&) and verify the 
resulting solution. The simplest example is obtained by assuming the 
incident wave to consist of unit displacement lasting for a finite time, say 
27, and putting, therefore, 


fg)=1 (él <7); f€)=0 (l€| > 7). (17) 
It then follows from (15) and (16), after some simple calculations, that 
the reflected wave is given by 


xsin Co ‘ , 
— oft—* sin y yen) 





, 


A : dog f=7 I = 
1—)? 2 
Toul ©) +a 


and the transmitted wave by 


n-af-t-3)} 


> 


g(€) 











(19) 
lean a T—€ 3TH Mog (P= Eta? 
ml Ta \°o" “—" +tan-! +5 S(T Le 7? 


It is easily verified that the boundary waits hold except when 








and these discontinuities can be disregarded.t 


+ The integral theorems expressing the conservation of momentum and the continuity of 
deformation continue to hold near the singularities. 
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The logarithmic term in (18), which characterizes total reflection, gives 
rise to logarithmic infinities corresponding to the discontinuous beginning 
and end of the incident pulse. It is an odd function of €, and represents 
an additional displacement gradually increasing and becoming infinite at 
the instant corresponding to the arrival of the reflection of the sudden 
onset of the incident pulse; it then swings in the opposite sense, becoming 
negatively infinite at the instant corresponding to the arrival of the 
reflection of the discontinuous end of the incident pulse; and finally dies 
down, approximately proportional to the time elapsed. 

The transmitted wave (19) is a disturbance propagated without change 
of shape along any plane y = constant with velocity c. The tan-! terms 
involving inverse tangents give an even function of € which can be looked 
upon as a distorted version of the incident pulse; the logarithmic term 
is an odd function of € which contributes an additional displacement 
consisting of two successive ‘waves’ in opposite senses. As the distance 
from the interface, y, increases, w, dies down as 1/y. This slow rate of 
dying down distinguishes the total reflection of a pulse from the case of 
an infinite harmonic wave train, for which the disturbance falls off 
exponentially with penetration depth. 

This solution of the problem of total reflection can also be applied to 
electromagnetic waves. Consider, for example, the reflection of a plane 
polarized wave at an interface between two non-conducting homogeneous 
dielectrics when the electric field is orthogonal to the plane of incidence. 
Taking this plane to be the ay-plane, this case arises when the only non- 
vanishing component of the electric field is E,, while the magnetic field 
H, = 0, Assuming that there is no dependence on z, Maxwell’s equations 
in the upper half-plane are 


ok, _Bé H, oE, pe oH, ae 








é ! j 0, 7 7 0, 
oy c at Ox c ot 

eH, eH, KéeE, 

Ox cy c at is 


where K, » are the specific inductive capacity and the magnetic permea- 
bility respectively, and c is the velocity of light in free space. We can 
then take 
. od c é ce 
E=—, ee. tn So, 
ot : pb Ox 


pe ey” 
provided that ¢ satisfies the wave equation 


o*h : aod a pk ad 


bat ae y > 0). 
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In y < Osimilar equations hold, with K’, »’ instead of K and wp. At y = 0, 
H,, and E, must be continuous; hence 

od c Od 

at’ pe Oy 
must be continuous. The problem is therefore formally identical with that 
which has been solved. 

The problem which has been considered is perhaps the simplest of this 
type. It is clear that a number of other problems involving plane waves 
of arbitrary shapes can be treated on the lines indicated. One example 
of some interest is the propagation of waves of arbitrary shape over the 
free surface of a semi-infinite elastic solid, analogous to Rayleigh waves. 

Assuming the motion to be two-dimensional, let the solid occupy the 
half-plane y > 0. The equations of motion of an elastic solid can be 
satisfied by taking a » ‘ 

; = op erp 7= Op Op (20) 
ax " dy’ ey ea’ 
where wu, v are the components of displacement in the x- and y-directions 
respectively, and ¢, % satisfy the wave equations 


Op, 0 186 YM A 12 


¢ ¢ oa 9 ayo? c 9 
dx? oy? —s a* Ot” da? * dy? BP 





oo" 
a 
a and 6 denoting the velocities of propagation of longitudinal and trans- 
verse waves. Now we can assume that 

d = f{x—ct, y(1—c?/a?)} = f(é,n), say, (21a) 

b = g{x—et, y(1—e* te = 9(€,7'), say, (21 5) 
provided that f(€,) and g(€, 7’) are plane harmonic functions. Then the 
displacement at the surface, y = 0, is propagated without change of form 
with velocity c. By giving this velocity an appropriate value (which must 
be less than b) the boundary conditions at the surface can be satisfied. 

In fact, the normal and tangential stresses which must vanish at y = 0 

can be written in the form 


— > Ov ou. av 
ae =— p{(at—208) +a? I, = phe(: +2), 


ex oy 





where p is the density. Using (20), (21a,b) and replacing é?f/én?, 2g/én” 
by —éf/a&*, —e*g/0&* respectively, it follows that the boundary conditions 
are 
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When » = 7 = 0 we can write 7 for 7’ and replace these equations by 
the equivalent conditions 


(s- 1) £46. 0)— (1-5), 0) = 0 (22.4) 


2b? 


(1 — 5) fal€. 0+ (Sp) oe 0) = 0, (225) 
a”, ' 


where partial derivatives are indicated by suffixes. 

We now make use of the fact that the derivatives h,(é, 0), h,(€,0) of a 
harmonic function are conjugate functions and constitute, if h fulfils 
suitable conditions, a pair of Hilbert transforms (2), satisfying 


he(é, 0) =. P h = dx. h,(€, 0) a Ip eos da. (23) 
° * gf ay 7 a—§ 


Formally, these relations can be derived by putting first 6 = 0 and then 
i = 0 in (14). Taking the Hilbert transform of (226) we obtain therefore 


c2\1 2 
(1 — 3) Se 0Q)— Case 7 9,(€, 9) = 0. 
a?j °° 2b? 


If this is to hold at the same time as (22a), c must be a root of the equation 


(s-)-(-3h0-3) a 
20 a y) 


This is the equation for the velocity of propagation of Rayleigh waves (3); 
it always has a positive root less than b. 
Let c now be the relevant root of (24). Then, at the surface, 


: é : o*\3 
Ug(F) u(e, 9) = fel€,9)+ ims 9,,(§, 9); 


> 


Vo(€) v(€,0) = (1S) a(€ 096.00 


By (22a), (226) these equations can be replaced by 


2 ‘ c2 
felé.0),  % = —T9¢lE,0). 


Uu 
0 3 me 
2b? 2b2 


bo] 6 


Hence, using (24), we find that the components of displacement at the 


iree surface are linked by the equation 


; (€ = x—ct) (25) 
oa x 


— C 


(2, 
5 E.. 4 P [ Ul) 7, 


tae 


and there is a corresponding equation for uy in terms of vy. It is apparent 
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that while these surface waves are propagated with the constant velocity ¢ 
and without distortion, uw) and v, cannot both have a sharp wave front; 
there must be some displacement at all times. 
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SUMMARY 


This paper discusses a solution of the laminar boundary-layer equations recently 
published by Professor 8. Goldstein. Conditions near an assumed singularity at a 
point of separation are investigated numerically, tables are given for general use, 
and a comparison is made with numerical information from other sources. In 
addition, an improved form of the ‘outer’ solution of Karman and Millikan is 
described. 


1. Introduction 

In a recent paper (1) Professor Goldstein has discussed an asymptotic 
form of solution in the neighbourhood of a separation point of a laminar 
boundary layer. The solution has some curious features connected with 
the determination of its arbitrary constants, and some doubt is, at first 
sight, thrown on its validity. Nevertheless, in an investigation which 
showed the need for a detailed examination of this problem, Hartree 
succeeded in obtaining very good agreement (within a range limited by 
the scope of his numerical results) between velocity profiles at separation 
given, on the one hand, by G{20]+ of Professor Goldstein’s paper, and 
on the other by his numerical solution for the whole boundary layer. 
In view of this agreement it seems desirable to coordinate the two investi- 
gations further, and to attempt to remove some of the uncertainty 
referred to above. 

Interest centres first on whether the integral condition G[35] is satisfied, 
and secondly, on whether the disposable constant a, can be determined 
by applying the suggested condition at separation. This paper seeks to 
examine these questions largely by computational methods, obtaining in 
the process tables suitable for applications. Whilst formal proofs have 
seemed too difficult to attempt with existing techniques, the results of this 
investigation contribute strong evidence that both questions may be 
answered in the affirmative, and that therefore Professor Goldstein’s 
solution is very likely to be valid. 

In examining the second of the two questions, it seemed desirable to 


} References in this style are to equations in Professor Goldstein’s paper. Definitions 
and notations, which are the same as in that paper, are not repeated here. 
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develop an improved version of the ‘outer’ solution due to Karman and 
Millikan, and a brief summary is given here of this extension; it is intended 
to discuss this more fully elsewhere. 


Numerical investigation of the solution 

The general plan of attack consists of computing in succession the 
functions fy, f,, fo, fs, fs, and f;—or, in the cases of f, and f;, those parts 
of the solution which are obtainable—and then using the results, firstly, to 
compute values of the integrand in G[35], secondly, to calculate the velocity 
distribution in the neighbourhood of separation, in order to make a direct 
comparison with Hartree’s values, and thirdly, to construct the velocity 
profile at separation and attempt to join it smoothly to the improved 
outer solution. 

It is preferable not to use the explicit forms for the functions f, deduced 
by Goldstein, except for the purpose of starting and checking the integra- 
tions. The method of numerical integration to be described uses the 
asymptotic expansions to start at some convenient large value of 7, and 
determines f;, f,, f,; by integration towards the origin » = 0. In this way, 
accumulation of error is minimized since there is no solution tending to 
become larger in the direction of integration, and rounding-off errors 
therefore die out rather than magnify themselves. Consideration of this 
is of some importance when each function is used in the determination of 
subsequent functions, and it appears in fact that, without any special 
precautions being taken, there is no loss in the number of significant 
figures as we proceed from one integration to the next. 

The method, which is becoming a standard one in other connexions 
(especially in the making of tables to many figures, where a powerful 
process with automatic checking is desirable), consists essentially of a 
step-by-step application of a Taylor series, the coefficients of which are 
recalculated by a recurrence formula at each point of the integration. 
If h is the length of a step, and @ a variable fraction, then a function f(7) 
may be expanded in the region of a — 1 = % in the series 


oe 
f(no+ 9h) = f(n9) + nf (No) + rl ()+-.-.-- 


For simplicity we writet 
— f(m9) = 7°f (no) 
and abbreviate the right-hand side to 7”. Then 


FS (no +Oh) = f (no) +07+ 6772+ 73+... . 


+ See British Association Mathematical Tables, part-volume B (The Airy Integral, by 
J.C. P. Miller), p. B7. 
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Also, by differentiating with respect to the independent variable 6h, 

Th (no +Oh) = hf'(notOh) = 7+207?+ 36773 +... , 

Tf (mo t+Oh) = $h?f"(not+Oh) = 7?+3073+ 66774+.... 
In applying these to the numerical integration we put 6 = +1, obtaining 

(moth) = f(n)Ar+ Pt A+ At 7°+... 

of (no--h) +r42724344744 54... }. (1.1) 
Tf (no +h) + 7?+373+ 674+ 1075+... 


If the function and its first two derivatives (the equations being third- 
order) are known at some point 7 = mp, then values for the next step can 
be found by using (2.1) in conjunction with the differential equation. 
The formulae (2.1) with the alternative minus signs are at the same time 
used as a check by reproducing the values of the previous step and so 
testing that no significant arithmetical error has arisen in the working 
of that step. For ease of calculation the terms of (2.1) are grouped as 
follows: 


f (mo +h) (79+ 7? +-74+-76+....)4-(7!+73 +75...) 
th (Noh) (71 4+- 373+ 575+...) +(27?+ 474+ 675+...) ’ (1.2) 


Tf I h) (7 


+- 674+ 157%+-...)+(373+ 107°+ 2177+...) 


9 


The quantities in brackets are calculated separately and added or sub- 
tracted according to whether the next step or the check is required. 

This method has numerous advantages in making tables, by comparison 
with standard finite difference methods, and even in general application 
in a problem such as the present one it can effect a great economy of 
effort. Its fundamental advantage is that it is, in principle, unlimited in 
accuracy when the length of step / is specified. This is not true of a 
difference method since high-order differences are inaccurate because of 
rounding-off error, and may in any case contain few figures; there need 
be no such loss of accuracy in the calculation of 7” for large n, hence in 
general longer steps may be taken and the computation appreciably 
shortened. Even in cases where complexity forces the method to be cur- 
tailed to include derivatives no higher than in the differential equation— 
and therefore to be equivalent to a difference method—the present writer 
has found it more convenient on account of its very direct checking 
procedure. 

The functions fy, f,, and f, are of course easy to tabulate directly. For 
convenience, in all that follows, the functions considered will be f,/a4 and 
will still be denoted by f,; this is effectively the same as putting a, = 1, 
or as replacing the other variable ¢ by a, é. 
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2.1. The function /,(7) 
The differential equation is 
3 —37°fs +305 —6nf,; = —10a, n?—$n° = —17-78480 3n*?—§y)?. 
(2.1.1) 
Translating this into reduced derivatives we have 
—7* = thy r?— fh? n?71+ h3n7r?+h3(2-96413 4n?+0-27)), 
with similar formulae, obtained by differentiation, for higher orders. On 


choosing a step h = 0-05 and rearranging into a form more suitable for 
computing, the first three become 
—7 = = a [20 {679+ 20(—3-57!+ 20y7?)}]+ 
960000 ; 
+0-00037 05167 y?-+0-00002 77), 


] 





_ 2 0 9 f e 9 , P) ‘ 3\) 
—7 = 67° + 20n} — 71+ 20n(— 47°+ 607 { 
384006 vi} T i\ + 607%) } | 
+-0-00000 92629 18y-+-0-00000 17361}, 
, ] 
<7? = — [2-53 + 20n{ — 57? + 20n(— 1-573-+ 120y7*)}]+ 
9600000. — 7( + 12074)}] 


-+-0-00000 00926 292+-0-00000 00694y3. (2.1.2) 
The asymptotic expansion of f, is 
3 ~ —§+7:11392 1n+3-31100 8y?+-0-95597 76?(4 log y— 4-60728 6)— 
‘ es - ~ Fes 15 1 105 1 
— $n*—0-25492 736] 4y°—— —4+— ——-—_ —_.1. ,,,], (2.1.3) 
8 y* 167% 32 71° 
The integration was started at » = +4-00, using the following values 
computed from (2.1.3) and (2.1.2), and checked by substitution back into 
(2.3.1): 
fs — 338-572 7 = —25-3090 7? = —0-°74695 
7r? = —0-01149 0 7? = —0-000100 75 = —0-00000 02 


It was taken by steps of 0-05 as far as 7 = 0:5, and checks were applied 
at convenient points using the asymptotic expansion. Finally, a second 
integration, starting at 7 = 0 and proceeding for a short distance in the 
positive direction, was carried as far as 7 = 0-70 with steps of 0-05, giving 
some overlap and check. The accumulated error after 65 steps of the 
main integration was found to be only 1 or 2 units in the last decimal 
carried for each derivative. 

The functions fy, /;, fo, fs are given, with their first derivatives, in Table |, 
and it is believed that the last figure quoted is correct. 
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TABLE 1 





‘ 7 o(7) f5(n) /1(m) }3() JTo(n) fo(n) } I3{n) I3{n) 

















0-0 0°00¢ 0°0 0000 0000 0°00 0°00 
(2.1.] I I o2 0-018 0°356 0°03 0°66 
2 I 004 «(04 0°071 o-7II O13 1°32 
3 0°045 06 0-160 1-064 0°30 1°97 
; ‘4 I 80 08 0284 1°414 0°53 2°61 

i 

' 

05 0-125 5 9 0°442 758 | 0-82 3°22 
rs. Or 6 180 I*2 0°635 2-091 I'l] 3°78 
ible for | 5 245 } I*4 0°760 2°410 1°57 4°22 

8 64 1°6 1-116 2-709 | 2°02 4°69 
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2.2. The function /,(7) 
The differential equation is 
fo —4nfa t+ 4n'fs — Ths = Py—60§ 9? — 2a 4° +4n°f3 — 16nfy + 12fy. 
(2.2.1) 
It is convenient at this point to define a new function 
Fala) = fal) 2141?’ 40 ab"), (2.2.2) 
since a, is unknown at this stage, and the effect of P, may be added in 
later for any given problem. Then the differential equation for f, may be 


written 


pu 


2 —h0fit+4n°fa—Tnf, = —18-97794n?—3-55696° + 47°fg — 16n fs +12/, 
= —A,(n). (2.2.3) 
The successive reduced derivatives are given by 


~ 960000 [20n{77° + 20y(— 471+ 20n7*)}]+ 3h3A4(n), 


i. 01 90+7f 


= sea000G"'* 20n{—7!+ 20(—5r?+-60n7°)}|+ th4A4(n), 


where the quantities on the extreme right are connected directly with the 
previous integration by the relations: 
$h?A4(n) = 0:00039 53737y?+0-00007 4103375+ 
+ (0-0627?—0-006y71-+ 0-00025r°)(—f,) 
xh*A4(n) = 0-00000 988437+0-00000 46315y4+ 
+-(0-05n?73— 0-0016n7?—0-00002 08371)(—f,). 
A,(y) is most easily computed from the differences of A,(n). The 
asymptotic expansion of f,(7) is as follows: 
fy ~ 0-06349 206y7—1-18565 35y4—5-3n3—17-11132 72+ 
+ 1529564 2y log n+-9-59998 567 — 
211,31 £3151 


1 
itt Hel. ~ te SS 
FE 7 2 27° 88 7? 


—3-82391 04 1 n+ 11 mm. = 35 a. . (2.2.4) 
5 ° 7* 4 "" 8 yt 
The integration was again started at +-4-00 using the following initial 
values: 
Jag = —1479-124, 7 — 105-6392, 7? = —3-2263, 73 — —0-05514, 


74 — —0-00057 8. 75 —Q-00000 3, 


and continued, as for /;, by steps of 0-05 down to y = 0-5. It was found 
that the contributions of f, and its derivatives to A,(y) were of sufficient 
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accuracy to give five or six significant figures in f, where there had been 
five or six in f;. This feature, which depends partly upon the correct 
choice of h appropriate to the accuracy required, persists in later stages, 
and means that the computations could, if desired, be repeated to a larger 
number of figures, with a smaller integration step, without any increase 
in complication. 


2.3. The function /,(7) 

The differential equation satisfied by f, is 

fs —4n°fs +3n°fs —8fs = Shi fa—#i Sat ht fat Sols —e fe + fe fo, 

(2.3.1) 

and, although the right-hand side of this is known completely (apart from 
the term involving a,), the complete solution for f; is not known. The 
equation (2.3.1) has three complementary functions »?, g;, h;, where, by 
G[52] and G[53], 





” =. (m —1)! (4)! l nimi 
Is (m-+-4)! (—42)! 4m—1 8m! (2.3.2) 
m= 
and h; 1+4n*—sts7°. 


Thus g, starts with a term in », and is the only complementary function 
containing exponentially large terms. h;, on the other hand, is a termi- 
nating series, and by its use it is quite simple to reduce the order of (2.3.1) 
and find that the required solution with a double zero is 


9 


fs(n) X50? + 4(a4+ ag Xg)(17—Gs5)—asFy 0° —sh9%3 99+ 








a 
ee | 1 bt" | H5(q)(—29+ §9>— aye" dy 
9 “~ , & 5 0 
TT 11 —34+--59-z ) —thesoeoapenperty aT dndn, 
J | ""s = 4—5*+ 837° —ts7 —17647" 
0 0 (2.3.3) 


where H,(7) is the part of the right-hand side of (2.3.1) which does not 
include the terms 
— 14(a4+- a a 3)9?—§P, 1° —§ag °9—HF, 0’. 


> 


The form in which (2.3.3) is presented warrants a little discussion. 
There is first the complementary function term a, 7? which in any case 
would have an undetermined coefficient at this stage. The second term 
associates the term in 7 of the particular integral with the non-terminating 
complementary function g,: the constant a, is at present undetermined. 
The fourth term has been introduced separately since it is the only other 
term involving «, explicitly, while the remainder, represented by the 
repeated integral, is a function whose Taylor series for small 7 commences 
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with a term in 7%. Use of (2.3.3) to determine f; would be prohibitive, 
as a brief attempt showed. On the other hand, the method of numerical 
integration used so far will not, without modification, give the full solution 
for f, since no complete asymptotic expansion is known. However, the 
following procedure may be adopted. 

The asymptotic expansion of the right-hand side of (2.3.1) is (apart 
from the term involving «,) known completely. This may be substituted 
and the equation solved in a series of descending powers of 7, starting with 
a term in 7°, to give an asymptotic representation of part of f;. This will 
not be complete since only one of the three boundary conditions is satis- 
fied; the asymptotic expansion so found ensures the absence of exponen- 
tially large terms, but does not necessarily correspond to a function with 
a double zero at the origin. It is therefore necessary to remove from it 
a constant and a multiple of 7, and to add an arbitrary multiple of 7? as 
usual. 

We may conveniently define a new function 


fs(0) = fs) —0%5 9? —4ag 9—A(1 +3 ahs 1) +P 1°, (2.3.4) 
where A is a constant to be determined later. When this, together with 
the asymptotic expansions (2.1.3) and (2.2.4) is substituted in (2.3.1), we 
find 

Fe’ —4n*fs +3n*fs —8nfs 
~ —1-90343 2n?7—6-12067 4n®—10-19709 475 log n+-7-5045775— 
—23-71307y!—57-1029973—95-21049n? log 7+ 246-381897?— 
—61-182577 log n—336-0986y—56-91137-+204-7121/+ 
+ 66-30731/n2—212-2270/n3— 163-7696/7°— 140-265 /7°+ 
+981-470/n7+.... (2.3.5) 
The solution of this, apart from an arbitrary 7? term (which may be taken 
as having been absorbed into a; 7”), is as follows: 
Fs ~ —0°47585 82n®—1-36014 98y5—2-54927 3n‘log n+ 2-51346y!— 
—9-4852373-+ 27-20300n log 7 —62-622547-+41-26757— 
— 10-23560/?—3-40037 5/3-+-5-89519/*—1-46223/7°+ 
+5-1006/n?—21-111/n8+.... (2.3.6) 
If an integration for f, is started and carried back to » = 0, the resulting 
function will have in its Taylor series for small » a constant and a multiple 
of ». The complete function f, is constructed by adding the four terms 
on the right of (2.3.4), hence we see that A is equal to the constant just 


mentioned, and that the coefficient of » must be equal to 4a,. Thus f; 
could be completely determined in this way (apart from the term a; 7) 
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were it not for the fact that the integration, by the time it reaches 7 = 0, 
has lost so many figures that the values of A and a, would be practically 
useless. However, it will appear that the contribution to the integral 
G 
relevant to the immediate argument. (2.3.6) was used to start an integra- 


35| of the terms involving A and a, is zero, hence this difficulty is not 





tion at 7 4-00 with the values 
f, 4269-28, + 277-087, 7? = —7-623, 7? = —0-1122, 
7 —0-00095. 
Only 7? needed calculation—7+* could be obtained from the differences of 


7—and was found from the formula 
l > 
K ¢ = ¢ . ok L¢ 2\i1_1 
7 = ——_ [20 n{81r°+ 20n(—4-5n7!+ 20n7?)}|]+ 3h3A,(n), 
960000 
where 


gh A;(n) (1579°7* —g0007* 4 subooT)( —f4) + 
+(isfe7°—atoSe 7 + sd00t2 7°) —fs)- 
The integration was continued as far as 7 = 0-5, and extrapolated to the 
origin by fitting a polynomial 
fs = a+bn+en*+dy>+O(n!) 
which corresponds with the known form of f, near the origin, and which 


gIVeS 


v= 81,, A = —0-84. 


The functions f, and f; and their first derivatives are shown in Table 2. 


3. The integral condition 

We now come to a discussion of the condition G[108], which expresses 
the fact that the particular integral of G[100] which has a double zero 
at the origin contains no exponentially large terms. For convenience, it 
may be restated here as follows: 


If the particular integral of 


fe —37°f6+5n"fe—9nfe = He(n) (3.1) 


satisfies the above conditions, then 


’ =e 1_.10)p—8 , 86 
H,()(» 57° +7607” )é id dy = 0. (3.2) 
0 
Evidently, isolated terms of H,(7) which demonstrably do not give rise 
to exponentially large terms in f, may be dropped from the right-hand 
side of (3.1) and hence from the integrand of (3.2) before formulating the 
condition, 
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The form of the condition raises an immediate difficulty; to prove such 
an expression equal to any specified number by computational methods 
is, strictly speaking, impossible. At most, upper and lower bounds can 
be calculated, but even this has a point only when it seems likely that the 
expression can be proved not equal to the given number. However, we 
may investigate its value within the limits of accuracy of the previous 
computations. If we write 


A=, fs =fatosn, fa = fetogn?+tA(n>—1857"), 
fs = Fotos 9° 40g +A 1+ $05 shan) —aP, 0° 
in the right-hand side of the equation for f,, this becomes 
— 16a; 9° — 16014(9 +a 9°+-349°) — 805 °— BP, 0° — for Pq? + 36") + 
+A(16— Ht — 888) + [40°F — 20nfs + 16f5+ Sfofs —10fs fa + 
7 Tfefa+ Offs + 1203 f+ bag nfs —5fs?— 20 x3 7f3 }. (3.3) 


The particular integrals corresponding to the first few terms are: 


4a, n+ 2 Jag n+ 4ar4(xg 7 — $4) — sal, 7° — derq Py n®+ $A(n?—shpn"), 


and to give this expression a double zero, it is sufficient to add on a suitable 
multiple of 96> the terminating complementary function, i.e. to add 
—(4a5+ 4ag «4+ 20§)(n+3470°—3357°). There is no constant to be removed 
by use of h, and therefore these terms induce no exponential terms in f,. 
The condition (3.2) may therefore be replaced by 


| Ie(n)(n?—308 + shone" dn = 0, (3.4) 

a 
where J,(7) is the quantity in square brackets in (3.3), which consists 
entirely of quantities which are known at this stage. The integrand of 
(3.4) may be calculated easily from the results of the numerical investiga- 
tion. It is a little unsatisfactory that there is some loss of accuracy in 
doing this since the various groups of terms tend to cancel partially within 
themselves; thus, for example, at » = 2 


4n?f; —20nfg +16f, = +9360—12990+2990 = —640. 


The final product of the numerical work is thus not as accurate as might 
have been hoped from the integrations, which were done using six or 
seven figures. The effective range of integration is from about » = 5 
to » = 3:5, above which the exponential term makes the integrand 


negligible. The form of the result is shown in Fig. 1, where the curve is 
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seen to have two positive and two negative regions. The areas determined 
by numerical integration are as follows: 


3°5 
| = —457-4 The fourth figure is believed to be reliable. 

2-340 

2340 
| = -. $87 Only three figures are known here. 

1561 

1561 
| =—1-7 The second figure is unreliable. 

1-425 

1-425 
| = +68 The second figure is definitely unreliable, but 
0 should not be more than 3 or 4 units in error. 

Sum = —4. 





| 


“1000+ 


eo 
Fic. 1. Values of the integrand of f J4(n) (4? —4y° +zh5y")e-t"* dn. 
0 


As regards the portion of the integral beyond 3-5, extending to infinity, 
the rapidly diminishing exponential makes this obviously small, but for 
a more formal statement we may consider the asymptotic form of the 
integrand. It can be shown that the part of the integrand apart from the 
exponential behaves for large » like 717; hence if a constant A is defined 


80 that ; 
Je(n)(q?—3°+ sh”) < An!” (all » > 3°5) 


the contribution to the infinite integral is less than 


A | net" dn = Ae-*89)2(3-5)14+4 56(3-5)!°-+ 1120(3-5)®+6720(3-5)?]+ 
_ +13440A,/(}7)[1—erf (3-52/2V2)] 
=0-7A (approximately). 
5092.4 D d 
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A brief examination of the computed values of the integrand shows that 
a coefficient A = 0-1 considerably over-estimates the quantity concerned, 
and hence the tail of the integral beyond » = 3-5 is quite negligible to the 
accuracy required. Thus it seems quite likely that the condition (3.4) is 
satisfied, and we may proceed to a more detailed examination of the 
velocity field. 


4. The velocity field at and near separation 


The results of §3 have contributed evidence in favour of there being 
a solution of the form considered. It appears likely, as Professor Goldstein 
has pointed out, that the solution contains just one disposable parameter 
—the constant o,—and that a value for this constant may be selected 
by setting up a smooth transition to the main-stream conditions at the 
section €= 0. If this is done numerically or graphically it does not 
follow, of course, that the value so found is uniquely determined by the 
condition, but it seems most natural that it should be. Hartree, in his 
investigation, estimated the value of «, by comparing the distribution of 
skin friction as given by Goldstein’s solution (as far as it was known at 
the time) with that arising from his numerical results. This indicated a 
value of about 0-43 for a,. A check was applied by computing the velocity 
profile at separation and again comparing it with the corresponding result 
from the numerical integration. However, as we shall see, this result of 
Hartree’s is in error. Also his comparison of separation profiles is of 
restricted value since it is confined to points quite near the wail, where 
the theoretical profile is remarkably insensitive to the value of a,. 

In this section we show how to derive tables for the velocity field at 
points near separation (the choice of «, being left open for the moment). A 
comparison is then possible between a computed velocity profile and one 
derived from Hartree’s integration for a certain point near separation. 
We then point out and correct the errors in Hartree’s calculations, and 
compare distributions of skin friction. Finally, we compare velocity pro- 
files for the separation point derived from different sources. 

The tangential velocity u, is given by 


oy 
U 


wu, = 
1 - 
oY, 


= 28 fot(a fi +(u €)*fe+---]; (4.1) 
where, as explained previously, it is convenient to remove the constant 
from each function f, and to combine it with €. (4.1) is in a convenient 
form for tabulation, and the quantity w,/2é*, derived from the contents 
of Tables 1 and 2, is shown as a function of 7 for several values of a in 
Table 3. The contributions of the terms involving P,, P,, etc., are 
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negligible, and the information in Table 3 is therefore effectively inde- 
pendent of the pressure distribution. However, it must be pointed out 
that Table 3 does not include details of the velocity distribution at the 
separation point itself, where the effect of the first term involving P, is 
of comparable magnitude with the contribution of f,, and where the 
variable » becomes infinite for all y. 

It is convenient here to compare the systems of coordinates used by 
Hartree and Goldstein. The former defines non-dimensional distances 
X, Y which are related to a representative velocity and length in the same 
way as in G[3|. The main-stream velocity at the leading edge is taken as 
a velocity of reference, and the representative length is chosen so that 
the main-stream velocity distribution is given by 

U = 1-}X, 
X being measured from the leading edge. Here, it must be noted, we 
have used capital letters X and Y although Hartree actually uses small 
letters. Then the two sets of variables (£, 7) and (X,Y) are related in the 


following way: 
X,—X\} 1/8—X,\hy 
é= Gor) 1=3(z—3) (4.2) 

On ah g 4ig 





and, from (4.1), 


u = $/{(X,—X)(8—X, fot (aH At(uasfet---], (4.3) 
where X, is the value of X at the separation point, determined by Hartree 
to be 0-9590 with an uncertainty of one or two units in the fourth decimal. 
The velocity wu in (4.3) is Hartree’s velocity and is related to Goldstein’s 
u, a8 follows: 

u = (1—}X,)u, = 0-8801 4. (4.4) 

It is evident that if a value of X be selected, at which the velocity 
profile is known from Hartree’s integration, then the corresponding pro- 
files given by Goldstein’s solution with various values of a, may be 
obtained quite easily from the information in Table 3. Interpolation to 
make the two solutions agree may give an estimate of the true value of 
1. The information given by Hartree in his report does not allow a very 
wide choice of points at which to try the fit. On the one hand, it is un- 
desirable to choose a value of X very near X, since the small uncertainty 
in the fourth decimal of the value 0-9590 is unnecessarily magnified, and 
a two-parameter fit would really be demanded. On the other hand, if X 
is chosen to be too small, ¢ will be too large for the series in (4.3) to be 
adequately convergent. Considerations of this nature led to the decision 
that X — 0-956 was the only suitable point; here € = 0-1437, 7 = 1-740Y, 
and u = 0-03633 ¥ (a, €)"f;. In Table 4 is shown the true velocity profile 
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at this section together with two profiles calculated from Goldstein’s 
solution and having a, = 0-418 and a, = 0-487 respectively. It is not 
obvious what is the best method of interpolating for «,, nor is it clear how 
far out from the wall it is wise to proceed. That is to say, for small values 
of Y the solution in series must be quite accurate but the numbers involved 
are small and the rounding off is significant. On the other hand, for larger 
values of Y there are more figures but the convergence of the series solution 
is in doubt. 
TABLE 4 


Velocity profiles upstream of separation 

















Values of u 
i True a, = 418 | a, = 487 
0-0 *0000 *0000 | *0000 
*OO15 “0014 “0016 
2 ‘0041 "0039 "0042 
3 "0079 "0075 “0080 
4 | ‘0127 *O122 0129 
os | 0186 0180 0188 
6 0256 0248 0258 
7 "0336 "0328 "0338 
8 0427 0418 0430 
9 | 0528 "0519 0532 
10 0639 0630 0644 
I 0761 0751 0765 
2 “0892 0883 0898 
3 | *1032 "1025 *1039 
4 1182 1176 ‘1189 
I'5 *1340 °1337 *1350 
6 *1507 | “1508 1518 
7 "1681 1687 +1694 
8 1864 1876 1879 
‘9 +2053 "2073 *2071 
X = 0956 
X, = 0°959 


It will be observed that for Y greater than about 1-5, the two profiles 
calculated from (4.3) do not straddle the true curve but both lie to one 
side of it. We may make a cautious estimate, then, that the series solution 
will not give four-decimal accuracy for values of Y greater than about 2. 


As regards the interpolation for «,, a little consideration suggests that 
Yi 

the integral fu dY (where Y, is some convenient limit) should be made 
0 

equal to the value derived from Hartree’s integration. If more complete 


data were available—i.e. for larger Y—the method would be equivalent 
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to making the displacement thickness equal in both cases. Taking ¥, = 1 
(not greater for the reasons outlined above) we find 
xX, = 0-47. 

There is a small point with regard to matching velocity profiles as we 
have done here. The velocity as given by the series (4.3) may be in error 
either because «, is incorrect, or because the convergence is too slow for 
the limited number of terms available. If, as Y increases, a marked 
sensitivity to a change in a, appears before the convergence is too slow, 
then some conclusion may be drawn as to what is the true value of a. 
If, on the other hand, the two effects are not easily distinguishable, as 
in the example under consideration, any estimate of a, is correspondingly 
less reliable. 

Nevertheless, the value a, = 0-43 found by Hartree seems inconsistent 
with the above examination. In fact, it turns out that his computations 
are wrong on this point, and may easily be corrected. From the formula 


r 


G[22] it follows that the velocity gradient at the wall is given by 
eu U, (eu wre —e ¥ 
voters => —) = 4U,{ 2a, €?+- arg EF + erg F4+-...] 
OY} y-9 —2N2\ey1/0 
00-8801 > eK Or 9 « <4 5 
————[ (xy €)2-+ 1-7785 (a €)3+3-3110(a, €)4+8-15(a, €)5+...] 
a 
(4.5) 
when numerical values are inserted. If we divide by £?, and return to 
Hartree’s variables for the right-hand side, we may write (4.5) in the form 
2(Cu/e Y), 


[100(X,—X) }}! 
+-0-12478{100(X 


s 


100 6-6337[ a, + 0-34526{100(X,—X)}4a?+ 
X)}taZ+-0-0596{100(X,—X)}4a#+...], (4.6) 
where, following Hartree, certain numerical factors are introduced for 
convenience. Comparison of (4.6) with the corresponding equation in 
Hartree’s report (2) reveals three errors; we have also an extra term which 
Hartree did not have. When values of (éu/éY), are computed from (4.5), 
plotted against €, and compared with the several isolated values found 
in Hartree’s integration, the result is as shown in Fig. 2. It seems now 
that this comparison gives a value slightly greater than 0-48 for a,, and 
the latter value is more in accord with the other estimate. 

We may note that Hartree, in making his estimate, plotted not (@u/@Y), 
but the quantity on the left-hand side of (4.6) which would become infinite 
at X,—X = O if the wrong value of X, were chosen. The result was that, 
whilst his curve was smooth for X < 0-950, a peculiar irregularity was 
introduced by the two points for X = 0-956 and X = 0-958 which did 
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not appear to lie on a smooth continuation of the main curve. There was 
not therefore a very good fit between the two types of curve. The irregu- 
larity is concealed in Fig. 2 since there we are not plotting an expression 
which is very sensitive to the trial value of X,. Hartree asserts that the 
‘bump’ is real, and gives some estimates of the possible errors to support 
this contention, but one would feel much more confident (a) if there were 
more points between 0-950 and separation, and (b) if the corrections near 
separation in the h?-extrapolation process were not quite so large. 


‘O7r e 
‘06 
‘OSE Ped 
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04+ 


“02 = 


‘01 a 











0 1 j 
0-1 E 0-2 0-3 
Fic. 2. Skin friction. 
The shape of the velocity profile near the wall at separation is given by 
U, = a, yi+a,yi+a,yi+.... (4.7) 


When we transfer to Hartree’s variables and substitute the known 
values of the a’s we find 


u = 0-05501Y¥?—0-00229 20a? Y4—0-00065 73a? Y5— 
—(0-00010 22af-+0-00000 477)¥&+0-00000 2903 Y7..., (4.8) 
or, dividing by the value of the main-stream velocity at separation, 


U ré . 2 , "A , 
T= 0-0625Y?—0-00260 42a? Y4—0-00074 69a3 Y5— 


—(0-00011 62a4+0-00000 543)¥&+0-00000 3303 Y7.... (4-9) 
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The profiles given by the last formula, with various values of «,, are shown 


bv broken lines 





in Fig. 3 for the region in which Y >2. They are 
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Fic. 3. Velocity profile at separation. 
Ore 
0-8F 
| 
e 
—— Improved outer solution 
aS 
© Howarth 
| 
Jj | ! —— 
3 4 Y 5 6 
Fic. 4. Velocity profile at separation. 


compared in thi 


s diagram with the separation profile given by Howarth (3), 


Table III, last column. The reliability of this column in Howarth’s table 


was not known, but the results collected together into Figs. 4 and 5 
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confirm its accuracy to the number of figures involved, and make it , 
convenient basis of comparison. In Fig. 4 it is compared with the im. 
proved outer solution of § 5 (see later), whilst in Fig. 5 the various calcula- 
tions for Y < 1-5 are plotted. In this region the right-hand side of (4.9) is 
effectively equivalent to its first term, to three-figure accuracy, and it 
matters little what the precise value of «, is. 











0-2r 
u / 
U —— Goldstein 
x Hartree 
© Howarth 
0-1 
ot 
a 
Pal 
rl 
il 
— ‘ ‘ aS a 
0 0-5 1 v r 


Fic. 5. Velocity profile at separation. 


Returning now to Fig. 3, it seems that a good fit is achieved when the 
value of «, is a little greater than 0-48; this is in good agreement with the 
revised estimate from extrapolation of the skin friction, and in fairly good 
agreement with the estimate from the velocity profile at X = 0-956, 
Professor Goldstein has suggested that the latter estimate (which gave 
x, = 0°47) is the most reliable of the three. 


5. An improved ‘outer’ solution 
In this final section we describe briefly an extension of the ‘outer’ 
solution of Karman and Millikan in the particular case of a linear pressure 
gradient (which has been used in testing all this work). We take the 
boundary layer equation in the form given on p. 165 of Modern Develop- 
ments in Fluid Dynamics, namely, 
oF on a= *('\- Ga) (5.1) 
a Ue U2] ey?’ 


where z = U?—vu?*, and the variable ¢ is defined to be 
} ue [ U(x) de. 
i) 


The notation here is that of Modern Developments, and the variables are 
the natural, dimensional ones. There seems no objection to this since 
reference to previous literature is made easier, and the connexion with 
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previous sections of this paper occurs at one point only. The left-hand 
side of (5.1) corresponds to the inertia terms, and the right-hand side to 
the viscosity term, in the ordinary equation with x and y as independent 
variables. The first approximation, taken by Karman and Millikan, is 
3 92 
Gz az - 
a7 = ve (5.2) 
Cp 
and this is equivalent to an equation in which the viscosity term is 
multiplied by U/u. A convenient second approximation is 
1 z\a% ” 
= y{ 1] —— —_}___., (5.3) 
2 U2) ap? 
and it is easily seen that this is equivalent to multiplying the viscosity 
term by 3(U/u+u/U). This quantity differs from unity by less than a 
given amount over a much larger range than does the ratio U/u, and the 


Nx 


Fs) 


a>} 
| 
-i 


solution stands to gain appreciably in accuracy. For the case under 
consideration 


U=b—b,x, ¢=bx—}b,2%, U2 = b8—2b,¢. 


The solution, to Karman and Millikan’s approximation, is 





2 = b%(1—erf £)4 2, | Ste (14 207\01—en)], (5.4) 
where (= a . (5.5) 
2 (vd) 
The form of (5.4) suggests solving (5.3) in the series 
bal fo(S) —Abfi(Z) +A*h2fo(C) —A*P*f3(2) + —...], (5.6) 
where A = 26,/b%. This gives the following family of differential equations: 
(2—fy) fo +40 fo = 9, (5.7) 


ou 


n—1 n—8s 
(2—fo)fn +40 fn—8"f, = > fe > (—1""-f, (n>1). = (5.8) 
s=0 r=0 


The boundary conditions are 
fol) =f,0) = 1; f,(0) =0 (r > 2), 
TA) == & (all r). (5.9) 
It is possible to simplify the numerical integrations considerably by a little 
analysis. Thus, if we write 2—f, = 2y in (5.7) we obtain 
yy’ +2ly’ = 0 (5.10) 
with boundary conditions y(0) = 4, y(oo) = 1. The solution of (5.10), if 
it exists, may be shown to be 


y = $[14+ Bv2 9,(fv2)+ (Bv2)*p9(fv2)+ (BV2)8p4(fv2)+...], 
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where 8 = y’(0) and 
9(£) = 4vrerfé, p(€) = }(1—e-*8’— va fe’ erf é), 
sé) = —hée-€(1+ fe-%*) + hvorerf (f+ e-%"— g%e-26*) — 
—/(da)erf(Ev3)+ ayeake(erf €)°(8—&). 
Thus by taking the limits as § > «0, 
Qi>4dv7, o2>4, 93> eV7—d/(47), 
we obtain an equation for 8 of which the first few terms are 
1/V2 = 0-88628-+ 0-353662-+-0-029763-+L ... 

and, with an estimate of y’(0) from this, the necessary trial and error 
consequent upon the non-linearity of the differential equation is very 
much reduced. 


The linear equations (5.8) may be simplified in another way. For 
example, when n = 1, 


(2—fy) fi +40fi—(8+ho)fi = —foifo, 


and writing 2—f, = 2y as before, 


Uf +2 f—(4—y")fi = 2y"(1—y). (5.11) 
It is found that the substitution 
h=9n 
effects a considerable simplification of (5.11), and indeed of the left-hand 
side of all subsequent equations in the series (5.8). We get 


ygi—2lg,—69, = —4€(1/y—]), 
in which, not only is the left-hand side simpler, but derivatives of y are 
no longer involved. 

In such a manner the equations may be prepared and integrated 
numerically. Then by substitution in (5.6) the velocity may be tabulated 
as a function of % for each value of ¢. A simple quadrature gives y, the 
distance from the wall, as 

y= zh 2vd {< (for each ¢). 


Nv Jou 

There is, however, one arbitrary constant involved here which implies 
that in plotting u as a function of y an arbitrary shift in the origin of y is 
permissible without, however, the shape of the profile being affected. 
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ON SOURCE AND VORTEX DISTRIBUTIONS IN THE 

LINEARIZED THEORY OF STEADY SUPERSONIC FLOW 

By A. ROBINSON (The College of Aeronautics, Cranfield) 
[Received 28 November 1947] 


SUMMARY 

The hyperbolic character of the differential equation satisfied by the velocity 
potential in linearized supersonic flow entails the presence of fractional infinities 
in the fundamental solutions of the equation. Difficulties arising from this fact can 
be overcome by the introduction of Hadamard’s ‘finite part of an infinite integral’, 
Together with the definition of certain counterparts of the familiar vector operators 
this leads to a natural development of the analogy between incompressible flow 
and linearized supersonic flow. In particular, formulae are derived for the field of 
flow due to an arbitrary distribution of supersonic sources and vortices. 

Applications to aerofoil theory, including the calculation of the downwash in the 
wake of an aerofoil, are given in a separate report (ref. 9). 


1. Introduction 
It is known that the elementary solution of Laplace’s equation in three 
dimensions—i.e. the velocity potential of a source in hydrodynamics and 
the potential of a gravitating particle in Newtonian potential theory—has 
a counterpart in the linearized theory of supersonic flow, viz. the velocity 
potential of the so-called ‘supersonic source’. However, the development 
of the analogy meets with obstacles which are largely due to the fact that 
the velocity potential of a supersonic source becomes infinite not only at 
the actual origin of the source but also everywhere on the rear half of the 
characteristic cone emanating from it. Thus, in trying to evaluate the 
flow across a surface surrounding a supersonic source, the resultant integral 
becomes infinite. This and other difficulties can be overcome by the intro- 
duction of the concept of the finite part of an infinite integral, which was 
first defined by Hadamard (1) in connexion with the solution of initial 
value problems for hyperbolic partial differential equations. A description 
of this concept is given in §2 below, and applications to problems con- 
cerning source and doublet distributions under steady supersonic condi- 
tions will be found in §§ 3 and 4. . 

The concept of a vortex in the linearized theory of supersonic flow was 
first considered by Schlichting (2), who obtained the field of flow corre- 
sponding to a ‘horseshoe vortex’ by a synthesis of doublets. Schlichting’s 
approach has been the subject of some criticism, as in certain respects the 
supersonic horseshoe vortex is different out of all recognition from its 
subsonic counterpart. However, it is shown in§ 5 below that more generally 
the field of flow due to an arbitrary vorticity distribution, under steady 
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supersonic conditions can be calculated in strict analogy with the method 
due to Stokes and Helmholtz in classical hydrodynamics. The results are 
in agreement with Schlichting’s for the particular case of a horseshoe 
vortex. 


2. The finite part of an infinite integral 

Let D(x,y,z) be an algebraic function of three variables so that the 
equation D(x, y,z) = 0 determines a surface = in three-dimensional space. 
The surface & divides space into disconnected components &,, in which 
D(x, y,z) is of constant sign; also, D(x,y,z) will be supposed to change 
sign across any ordinary point of &. Further, let f(x,y,z) be a real func- 
tion defined in a certain region R such that 

fix.y.2) = g(x, y,2z)+go(x, y, z)D-™+-9, (x, y, z)D-i" +14... + 

T Ji, Y; z)D- ++ Ik +1(%, Y, Z)D**+-...+9pim(2; Y, z)D=-..., (1) 

where n is a positive odd integer, n = 2k+1, k = 0, 1, 2,..., and the func- 
tions g(x, ¥,Z), Jo(@, y, 2),... are either all analytic everywhere except on &, 
or at least have derivatives of a sufficiently high order, which are bounded 
in the neighbourhood of =. At any rate it is assumed that the analytic 
expressions for these functions may be different for the different ~,,. 

Given a small positive quantity e, we denote by N(e) the set of all points 
s which satisfy an inequality |s—s,)| < e for at least one point s, of =. 
We denote the boundary of NV(e) by B(e), and we denote by R(e) the region 
obtained by excluding from RF all the points of N(e). Furthermore, given 
a curve C, a surface S, or a volume V in R, we denote by C(e), S(e), and 
V(e) the subsets of C, S, and V respectively which are obtained by the 
exclusion of the points of \V(e). 

The concept of the finite part *J of a (finite or infinite) integral J— 
where J is any line, surface, or volume integral of f on C, S, or V respec- 
tively, e.g. 


. 


(fdx, [fdedy, | fdedydz, 
C § p 
where C, S, and V are supposed to be bounded—will now be defined as 
follows. 

Given the formal expression for J, we denote by J(e) the corresponding 
integral taken over C(e), S(e), or V(e) only. Subject to the specified 
conditions of regularity, J(¢) will be finite and of the form 

J(e) = a,e*ti+...ta,_,e*+0(1), (2) 
where O(1), as usual, denotes a function which remains finite as « tends 
to 0. We then define *J by 


*J = lim (J(e)—aye-i"+1—...—a,_, 7). (3) 
e—0 
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As stated in the introduction, the concept of the finite part of an infinite 
integral is due to Hadamard (1), whose definition, however, applies to 
a more restricted type of integral only. Hadamard writes [J instead of *J 
which is used by Courant and Hilbert (3). 

It will be seen that if J is finite, then *J = J. Also, the finite part of 
an (infinite) integral is invariant with respect to a transformation of 
coordinates, provided the Jacobian of the transformation does not vanish 
on 2. In particular, if we are dealing with the finite parts of integrals 
involving vector quantities, the result is independent of a rotation of 
coordinates. 

There will be no occasion for confusion if in future we refer to the finite 
part of a (finite or infinite) integral simply as ‘a finite part’. 

The finite parts of m-fold integrals in n-dimensional space, n > 3, 
m <n, can be defined in a strictly analogous manner. The rules valid 
for them are, mutatis mutandis, the same as for finite parts in three 
dimensions. 

The rules of calculation with finite parts, such as the rules of addition, 
are the same as for ordinary integrals. Also, if f depends on a parameter A, 
but D is fixed, then—provided the g functions are sufficiently regular (e.g. 
if they are analytic in the various &,,)—it is not difficult to show that we 
may differentiate under the sign of the integral, e.g. 


* * 
d fe. 
a [sae] = f Za (4) 
c c 


Under similar conditions the finite part of a multiple integral may be 
obtained by successive integration (including the operation of taking the 
finite part) with respect to the independent variables involved, taken in 
any arbitrary order. Thus, with the appropriate limits we have, for 


instance, 
* 


i f dxdydz = | ( i j faz) iy) dz. (5) 


More generally, we shall encounter cases where D, and therefore %, 
depends algebraically on one or more parameters. We shall show (i) that 
even in that case we may ‘differentiate under the integral sign’, and 
(ii) that if a given integral, or finite part, involves integration with respect 
to such parameters, as well as with respect to one or more of the space 
coordinates, we may exchange the order of integration without affecting 
the value of the integral. 

To see this we increase the ordinary three dimensions of space 2, y, 2 
by the parameter or parameters involved. Then, in the augmented space, 
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the surface D = 0 is again fixed, and in order to prove our assertions, it 
is sufficient to show (i)’ that in order to find the derivative of a finite 
part in n-dimensional space with respect to any variable which is not 
involved in the integration, we may differentiate under the sign of the 
integral, and (ii)’ that for any multiple integral in n-dimensional space, 
l1<m <n, we have 


.) = 
| (| fdry) dey..de — = | faeydey. AX», 


taken over the appropriate regions. It is clear that (ii)’ will prove (ii), by 
induction. 
We may reduce (i)’ to (ii)’. In fact, (i)’ states explicitly that 


* * 

Fr, of 
—- f dx, - ALn_s = [ — dx, ...dXy_1, 
O24 m @ ny Ox m 





and this will be proved if it can be shown that 


* 
. 


* 
f f daz,..dz._, = [ (| As dz,...d2,,. 5) aC 


where the lower limit of the integral with respect to z,, is arbitrary and C 
is independent of x,,. Putting 


. é 
F=— f . 
OX 
: 
we have f=fot F dz,, 


where f, is the value of f for an arbitrary but definite value of z,, (for 
given 2,,...,%,_—), and the integral is taken with that particular value of z,, 
as lower limit. Now, assuming that (ii)’ has been proved, we have 


* , * * , * 


[ | [ F diy) dey. dx a-i = { ( [ F diy dy. dx, 


* [ * 


and so (f—fo) da,...dx,, {Sz ——d2,...dX 1] AXm, 


. 





om” 
f dz,..dz,_,= f | f, Fa, AX, 4 | AX) + [ fod. 


and the last term is independent of z,,, as required. 
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To establish (ii)’, we have to prove that 


* * * 
| ( { fay) dz,..df_, = | fdey...dey. (6) 
* 
Putting | fda, = F, we see that (6) becomes 


* 


* 
| F Ax,...dLm1 = | SO by. (7) 
5 R ” 


taken over a certain region R on the right-hand side and over its boundary 
S on the left-hand side respectively. This is essentially the theorem of 
Gauss (or Green) for higher spaces. For m < 3, this theorem will be proved 
below for finite parts (without relying on the results of the present 
discussion), and the proof for greater m is quite similar. 
An important example of a finite part will now be calculated. Let 
D(x, y,z) be defined by D = a?—f*(y?+-2?) and f(z, y,z) by 
ox 
f(x,y, 2) [a — Pye pe) 


and 


for x? > B(y?+27), «x >0 





f (x,y,z) = 0 elsewhere, 

where o and f are arbitrary constants. We find that all the conditions 
laid down at the beginning of this paragraph are satisfied in every region 
R not including the origin, the surface = being given by the cone 
a? — B(y?+-2?) = 0. 

Further, let the open surface S be given by x = a, y?+2? < r*, where 
a >0O and r>a/f. S is a circular area including the circle x =a, 
y?+2? = a?/B8?, on which f becomes infinite of order 3. We are going to 


4 
evaluate *J = | f dydz. 


S 
Given « > 0, let S,(e) be the points of S(e) for which a? > B(y?+2’) 
and S,(e) the complementary set of S(e). Then f vanishes on S,(¢) and so 


Ile) = | fdyde= [ fdydet+ | fdyde= | fdyde 








S(e) Sie) S:(e) Si(€) 
27 pn/a 
x oa dydz _o p dpdé 
~ Ry BY J 
Si(e) 0 0 


where y = (a/f)pcos@, z = (a/B)psin@, and 7 is the radius of the circle 
bounding S,(e). It is easy to deduce from the definition of S(e) that 
(compare Fig. 1) 


1 = gla—e(1 +B}. 
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Ile) = | ] |: “1  Qaro 1 1 
) a gee | ae =a eS ] 717 
Plo, "Fp ey 
a — 
= al v{(2,)(1-+B%)}-He4-14 O(e4)] 
o) ‘ j 
1] ox dydz 2n0 
and so ‘J = | “heer ewer 7 eee’ ee (8) 
J [2°—By?+2*)}F p? 
This result is of fundamental importance for subsequent developments. 


Az 














Fia. 1. 


We shall also require extensions of the divergence theorem of Gauss 
and Green and of the curl theorem of Stokes to finite parts. Particular 
cases of the divergence theorem are in fact proved and applied by 
Hadamard in the above-mentioned book. 

A function f will be called an admissible function if it satisfies the 
conditions laid down at the beginning of this section. A vector function 
f will be said to be admissible if its components are admissible—and this 
is independent of the system of coordinates. If all the first derivatives 
of the components of f are admissible, then divf also is admissible. We 
are going to show that under these conditions we have for any volume V 
bounded by a surface S such as is considered in the ordinary divergence 
theorem, 


(9) 


* * 
| f.dS = | divtar. 
iS V 


5092.4 


Ee 
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By equation (1), the vector f = (f;, fo, f;) can be divided into two parts, 

F = (f,, f,, F;), and G = (G,, G,, G;), f = F+G, so that the F, are finite 

and continuous on = while the G; become infinite there. Then divF js 

either finite and continuous everywhere in its domain of definition or it 

becomes infinite of order 4} on &. Even then | div F dV exists as an 
V 


ordinary improper integral, and J F.dS = [ divF dV. Since 
§ Vv 


divf = div F+div G, 
it is therefore sufficient to show that 
* * 


| Ci = | div G dV, 
BA a V 
i.e. [ G.ds— { divGav = 0. (10) 
s V 


Let S’(e) be the join of (the set of points common to) V and B(e). Then 
V(e) is bounded by S(e)+S’(e), the union of S(e) and S’(e). Hence, apply- 
ing the divergence theorem to the volume V(e), we obtain 


| G.dS+ | C.d8 = | divGav 


ste) S'(6) Vie 
or | G.dS— [ divGdV = — | G.dS. (11) 
Sle) Vie) Se) 
Now [ G.dS is of the form 
Sie) 
* 
[ G.dS = aye“ 4. pay ett | G.dS+H(6), 
S(e) S 


where H(e) is a function which tends to 0 as ¢ tends to 0, and similarly 
f div G dV will be seen to be of the form 
V(e) 
* 
| div G dV = b,e-#"+... +b, e#+ | div G dV +K(e), 
Vie) v 
where K(e) is a function which tends to 0 as e tends to 0. 


* 
In other words, { G.dS and | divG dV differ from | G.ds and 
Sie) Vie) 8 


* 

| div G dV respectively only by vanishing functions of ¢ and by fractional 
; 

infinities of «. Hence, in order to prove (10), it is, by (11), sufficient to 
show that 


I G.dS = cye-"+... +0, ++ L(€), (12) 


Se) 





whe 
bet 
tion 
of ti 


whe 
holc 

N 
let . 
ord 


and 


wa) 


whi 


wh 
sid 
fini 
twe 
suf 


parts, 
finite 
vF is 
or it 


as an 


(10 


Then 
ipply- 


(11 


vilarly 


5 and 


tional 


ent to 


(12) 











LINEARIZED THEORY OF STEADY SUPERSONIC FLOW 415 


where L(e) is a function which tends to 0 as e tends to 0. And (12) can 
be readily deduced from the fact that the components of F satisfy condi- 
tions of the type indicated by (1). In fact (1) shows that on S’(e) G, is 
of the type 

G, = Que" +...+C, 74, (13) 
where the C depend on the parameters of S’(e), and similar expressions 
hold for the other components of G. 

Next, let f be a vector function of the same description as before, and 
let J be an open surface bounded by a curve C such as is considered in the 
ordinary curl theorem. Under these conditions we are going to show that 

? 4 
| f.dl = | curlf.dS. (14) 
C S 
Splitting f into two parts F and G as before, we first show that 
* 


* 
| G.dl = | curl G.dS. (15) 
C S 


Let C’(e) be the join of S and B(e). Then S(e) is bounded by C(e)+C’(e) 


and so, applying Stokes’s theorem to S(e), we obtain 


G.dl+ { G.di= [ curlG.ds. (16) 
ole) o%(e) Sle) 


In order to be able to deduce (15) from (16) we have to show, in the same 


e 


way as in the proof of the divergence theorem, that 
[ G.dl = cge-"+...4c,-44 Le), (17) 


C’(e) 


where lim L(e) = 0, and this follows from (1), as before. 
e—0 


We still have to prove that 
* * 


( F.dl = [ curl F.dS. (18) 
C S 


This is obvious, by Stokes’s theorem, if curlF remains finite every- 
where, and if curl F becomes infinite on = (in which case the right-hand 
side of (18) is an ordinary improper integral), provided S has not got a 
finite area in common with ©. Assuming on the contrary that S has a 
two-dimensional subset S bounded by C in common with &, it is then 


sufficient to show that 
on a 


( F.dl = [ curlF.dS. (19) 

C S 
Again, since curl F is an admissible vector, it follows that F is of the 
form F = F442), where the components of F® are finite and continuous 
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on 2 and the components of F® vanish on & (so that |F®|e-? remains 
bounded as ¢ tends to 0). Hence 


* 
[ F.al= | Fd. 
C C 
On the other hand, by the definition of the finite part, 
* * 
[ curlF.dS = [ curl F_.dS, 
SS rs 
and by Stokes’ theorem for finite functions f F® dl = J curl F®.dS and 
: . CG s 
so [ F.dl = { curlF .dS, as required. 
C s 


Equation (14) is now established completely. 


3. First applications to the linearized theory of steady supersonic 
flow 
In this and the following section we are going to discuss solutions of 
the equation 


a9, fay) a2, 
eo e®D aD 
ae a 
oz? 








—B = 0 (20) 
in relation to the linearized theory of supersonic flow. The following 
details are not intended as an exhaustive introduction to that theory; 
their purpose merely is to establish and explain the terminology used in 
the sequel. 

Assume that the free stream velocity of the given field of flow is parallel 
to the positive direction of the z-axis and is of magnitude U, where U is 


Ga? © dy? 


greater than the speed of sound 
- |dp 
@¢= =. 
dp 
Let the total velocity components in the direction of the 2-, y-, and z-axes 
be uw, v, and w respectively; then, assuming that u is large compared with 
v and w, and compared with its difference from U, we obtain, for steady 
conditions, the linearized Eulerian equations 
1 Op , ou 
BE PO poten 


E c ? 
p 0x Ox 


_l@ _ ya 





= : (21) 
p cy Ox 
lop ,,ow 
paz  ~ ex’ 


where the terms of second-order magnitude have been neglected. 
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Under the same assumptions, the equation of continuity which, in full, is 


‘ ‘ ‘ ' 9 
Gx Oy oz a? 


Ou. ov. dw. lfudp . vép . wep 
o,f 00, 1/0 ef ei) ay 
p Cx p cy p 


becomes, taking into account (21), 


( Ov . ow 
etd le ise (23) 


f ‘ ¢ , a 

Ox Oy oz 
where [? M?2—1] U2/a2@—1, M = U/a being the Mach number. 
Equation (23) is the linearized equation of continuity. If, in addition, 
the flow is irrotational, then we have curl q = 0, where q = (u,v, w), ie. 


ow ov ou ow ov Ou 
ose ean eae ane a Foy (24) 
OY Oz Oz OX Ox oy 
In that case, there exists a velocity potential ® so that q = —grad®. 


Expressing uw, v, w in terms of ® in (23), we obtain (20). 

Equation (22) expresses the fact that divq’ = 0, where the ‘current 
vector’ q’ is defined by q’ = pq. Since (23), which is the linearized version 
of (22), indicates that div|p(—f*u,v,w)] = 0, it will be seen that the 
corresponding ‘linearized current vector’ is q’ = p(—f*u,v,w), where p is 
now constant. Dividing q’ by p we obtain a vector q* = (—f?u, v, w) 
which will be called the reduced current velocity, or, briefly, the c-velocity 
of the flow. Thus, apart from the flow of q across a surface S, { q.dS, 
we are led to consider also the flow of q’ across S. In order to distinguish 
between the two types of flow, the flow of q’ will be called ‘c-flow’. It will 
be seen that the linearized equation of continuity (23) is the differential 
expression of the fact that the total c-flow across a closed surface 
vanishes. 

It now becomes natural to introduce alongside the conventional opera- 
tors V, grad, div, and A, the operators Vhf, gradhf, divhf, AhB. (Read 

hyperbolic nabla of index f’, ‘hyperbolic gradient of index f’, etc. The 
index 8 will normally be fixed throughout and may therefore be omitted.) 
The operator Vhf will be defined by 


me (- 2.2.3) 
Cx Cy Oz 


and gradh8 and divh§ as the two modes of this operator which apply 
to scalars, and to vectors in scalar multiplication respectively. The 
operator Ah £ is then defined by Ahf8 = div gradhB = divhfgrad. Equa- 
tion (20) may now be written 


Ah® = div gradh® = divhgrad® = 0. (25) 


\ 
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By the divergence theorem we have, for any function which is suffi- 
ciently regular on and inside a closed surface S bounding a volume J, 


[ gradh®.dS = f Ah@ar. (26) 
8 V 
If, in addition, ® is a solution of (25) then 


[ gradh®.dS =z @. (27) 
& 

However, if we replace ordinary integrals by finite parts, then equation 
(26) holds even when infinities are involved provided gradh® and Ah® 
are admissible functions with respect to a certain surface &. Thus, in that 
case, 


j gradh®.dS = -f Aho av, (28) 


while (27) becomes f gradh®.ds == @, (29) 


In the sequel, such surfaces = will be frequently composed of cones 
of the form («—29)?—f*[(y—yo)? + (z—Z»)?] = 0. 

We notice that if a function ¥ is admissible with respect to a surface 2 
and another function ® is regular (or has derivatives of sufficiently high 
order) on &, then the product ‘Y® is admissible. (The product of two 
admissible functions is not in general admissible.) We also notice that 
if a finite number of functions are involved, admissible with respect to 
different surfaces =, then they will all be admissible with respect to one 
and the same surface, viz. the sum of the different &™, given by the 
product of the functions D™ defining the &™. 

Let ® be a function which is regular (or has derivatives of a sufficiently 
high order) on and inside a volume V bounded by a surface S, and ¥ a 
function which, together with its first and second derivatives, is admissible 
in V (with respect to some algebraic surface X). Applying (9) to 
f = ¥ gradh®, we obtain 

* * * 
| ‘¥ gradh®.dS = ( grad ¥ gradh® dV + [ ¥Ah® dV, (30) 
S v r 
and similarly 
* * * 
| ®gradh ¥ dS = [ grad® gradh'¥ dV+ { OAh'Y dV. (31) 
S 


V V 


Now 
grad ¥ gradh® = —f? — 
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Ox oat ey Oy cz & 


= gradh ¥ grad®. 
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Hence, subtracting (31) from (30), we obtain 


* 


* 
| (¥Ygradh O—® gradh¥).dS = | (YAh ®—®AhY) dV. (32) 
S V 


This is the counterpart of Green’s formula, extended, however, to 
include finite parts (compare refs. 1 and 3). 


4, Source and doublet distributions in steady supersonic flow 
Elementary solutions of equations (20) or (25) are the functions® ,(z, y, z) 
defined by 

C 
J [(w—a9)*—B*{(y—Yo)* + (2—2)}] 
for (x—x,)* > B'{(y—yo)?+(z—z2)?], x > ao 


and ,(x,y,z) = 0 elsewhere, 


®,(x,y,2z) = 
(33) 


where P = (a, Yo, 2) and o are arbitrary. ®, will be said to be the velocity 
potential of a source of strength o located (or, ‘with origin’) at P. The 
actual velocity potential of a source travelling with a velocity —U, in 
sframe of reference travelling with the source, is obtained by adding — Ux 
to®, as given by (33). For reasons of simplicity, ®, as in (33) will be 
called a source for positive as well as for negative o. 

Similarly, a function ‘Y, will be said to be the velocity potential of a 
counter-source of strength o located at P if it is given by 

Co 

v[(e7—2%9)?—B{(y—Yo)? + (2—25)?}] 

for (x— Xo)? > B*[(y—Yo)? + (2—Z)?], L< Xp 


nd VR(z,y,z) = 0 elsewhere. 


¥p(x, y, Z) 
(34) 


A ‘doublet’ is obtained, by definition, by differentiating ®, with respect 
to length in any given direction, the differentiation being performed 
lative to the coordinates of P. Thus the velocity potential ©, of a 
loublet whose ‘axis’ is parallel to the z-axis is given by 

[(x—ay)? B*{(y—Yo)? + (z—Zp)*}]! 
for («—2y)* > B°[(y—yo)?+(2—2)?], & > Xo 


and O,(x,y,z) = 0 elsewhere. 


(35) 


A ‘counter-doublet’ is obtained, by definition, by applying a similar 
eration to ‘’,. An asterisk will be employed to indicate fundamental 
‘lutions of unit strength (e.g. 0%). 
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It will be seen that the velocity potentials of sources, counter-sources, 
doublets, etc., and all their derivatives are admissible functions in al] 
regions excluding their ego origins, the surface & being given by 

(w—2y)?—B*{ (y—yo)? + (z—Z)?] = 0. 
Also, the potentials of sources and counter-sources tend to zero of order 
4 as the affix tends to infinity in any direction not asymptotic to =, and 
seaitealy doublets and counter-doublets tend to zero of order $ under the 
same conditions. 

The velocity potentials due to line, surface, and volume distributions 
in — outside the distributions are obtained by evaluating the integrals 
| o® dl, [ of dS, | o@% dV, where o denotes the (variable) line, surface, or 
Vv eth density, and ®% denotes the velocity potential of a source of unit 
strength the coordinates of whose origin coincide with the variable(s) of 
integration. For sufficiently regular distribution functions o (e.g. if o has 
continuous bounded first derivatives) these integrals exist as ordinary 
improper integrals. For instance, for a surface distribution we obtain 


e a 
< \[(w--2%o)? 2 Br (y (y—Yo)* +(z—2»)?}]’ 


where o is defined as a function of the parameters u, v, of the surface 8, 


O(x,y,2z) = (36) 


given by 2% = 2% (U,V), Yo = Yo(U, VL), % = 2%(u, v) and the integral is taken 
over those parts of the surface for which 
, )2 2 . —— 
(w—ay)® > B*[(y—Yo)?+(z—2)?] and a > ap. 

The position is different in the case of line, surface, or volume distribu- 
tions of doublets, since the integrals corresponding to such distributions, 
viz. | oF dl, [ o} dS, and | oP} dV, are in general infinite. Thus, for 
a surface distribution of doublets whose axes are all parallel to the z-axis 

Y 
[ ; ae —oa(z 2—2»)B* dS = (37) 
2f 2 22h }8 
[(w—a'9)?—B*{(y—yo)? + (2—2)?} } 


which is, in general, infinite. However, provided c is sufficiently regular, 


we obtain 


the finite part of the integral still exists, and we may say that 


* 
© = [ o*dS 


is the potential due to a doublet distribution over S (with similar definitions 
for potentials due to line or volume distributions). An alternative method 
which avoids the use of the finite part and which has been used by 
Schlichting (2) and others, is to consider first the corresponding integral 
for sources (eqn. (36)) and then to differentiate with respect to (—?) 
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From a physical point of view this means that we calculate the potential 
due to two infinitely near source distributions of opposite strength. The 
final result is the same since, according to the rules given in the preceding 
section, finite parts can always be differentiated under the sign of the 
integral. It is precisely the possibility of carrying out all the necessary 
operations directly, without fear of obtaining meaningless symbols, 
which makes the finite part such a convenient concept. It will be seen 
that the alternative method is applicable only when all the doublets have 
parallel axes. 

In order to define the potential due to a volume distribution of sources 
(or of counter-sources) we have, as in classical theory, to take recourse 
to a limit process, as the integrand o@, tends to oo of the order 1 on 
approaching the point for which the potential is calculated. We therefore 
surround the point by a small sphere of radius ¢, evaluate the integral 
excluding the interior of the sphere, and then let e tend to 0. For finite e, 
the integrals in question exist as ordinary improper integrals, and the 
limit exists, as « tends to 0, since the volume of the sphere tends 
to 0 as e°. 

We are now going to show that the c-flow—defined as a finite part— 
across a closed surface surrounding a source of strength o is equal to 27a. 

We have to prove 

* 


* J gradh®,.dS = 2z, (38) 


8 
> 


where ®, is defined by (33). 


We may simplify the problem without loss of generality by assuming 


Ly = Yo = % = 0. Now let S be a small cylindrical surface bounded by 
two planes x +a and by the cylinder y?+z? = r?, where r > a/8. Then 
the integrand of (35) vanishes everywhere on S except in the circular area 
belonging to the plane x x. Hence *J reduces to 

* 


*J [ i oop dyda, ; 
J [a?—f*(y?+z*) |! 
Y*+2*<re 
and therefore *J 270, by (8). This confirms the theorem for the 


particular case of a circular cylinder. 

Next, let S be an arbitrary surface surrounding the source, then we 
may find a small cylindrical surface S’ of the above description inside S 
and we only have to show that 
* * 
[ gradh®,.dS = { gradh®,.dS. 
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Let V be the volume bounded by S and 8S’. Then by the divergence 
theorem for finite parts, (9), 


* * * 
| gradh®,.dS— | gradh®.dS = | div gradh®, dV. 
S S V 
* * 
But | divgradh®, dV = [ Ah®,dV = 0, 
v v 
since ® satisfies (25) and V does not include the origin. This proves that 
(38) is true generally. 
Similarly, we obtain for counter-sources, whose potential ‘Vp is given 
by (34), * 
— { gradh ¥p.dS = 2ne. (39) 
Ss 
More generally, if a surface S surrounds a finite number of sources of 
strengths o,, superimposed on an arbitrary field of flow which is regular 


inside S, then 
* 


— [ gradh®.dS = 27 }oa,. (40) 
S 

There is a similar theorem for counter-sources. 

In fact, (40) follows immediately from (29) and (38). 

We may also deduce, taking into account (i) early in § 2, that the c-flow 
across a surface surrounding a doublet vanishes, and more generally that 
the flow across a closed surface is not affected by the superposition of 
doublets either inside it or outside. 

Finally, it follows from (ii) early in § 2, that (38) can also be applied to a 
continuous distribution of sources inside a surface S, so that 


* 

~ { gradh®.dS = 2a [ 9, (41) 
S 

where | o is the total strength of the sources enclosed by S (and given as 
a line, surface, or volume distribution). The same applies in the limit when 
the distribution is actually bounded in parts by S. 

Equation (41) shows that the finite part of an infinite integral is more 
than an artificial analytical concept, and that in certain cases it may have 
a definite physical meaning. In fact, the product of density and c-flow, 
—p i} gradh®.d§, is the linearized expression for the total flow of matter 
whenever that expression exists, for instance if ® is given by a homo- 
geneous volume distribution of sources over the interior of a small sphere 
of radius « and centre P inside S, together with — Uz corresponding to-the 
free stream velocity. If o is the total strength of the distribution, then 
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according to (41), —p [ gradh®.dS = 2mop. Thus, 2zoap is the rate at 
which matter is produced inside S. Now let ¢« tend to 0, while o is kept 
constant. Then ® tends to —Ux+®p, where ®, is the potential of a source 
of strength o located at P. But o having been kept constant, it follows 
that the rate at which matter is produced inside S, and hence the rate at 
which matter crosses S, is still 27op. And this, by (38), can be expressed by 


* * 
—p | gradh®,.dS —p | gradh(—Ux+,).dS, 


so that the finite part is the natural generalization of an ordinary integral 
when the latter diverges. 

Applying (41) to a small surface surrounding a point P inside a volume 
distribution of sources, we obtain 


gradh® .dS = 2n lo dV, 


* 
8 V 
and transforming the left-hand side by means of the divergence theorem, 


this becomes 
* 


— { div gradh® dV = 2a [ odV. 
4 Vv 

Since this is true for an arbitrary small volume containing P, we must 

have 
Ah® = div gradh® = —2z0, (42) 
which is the counterpart of Poisson’s theorem in subsonic theory. 

Conversely, given the differential equation (42) over a certain region R, 
a particular solution of it is 

© = [ obfav. (43) 
R 
The general solution, as is easily seen by subtraction, then is 
© = [{ o*dV+4, 
k 
where ® is an arbitrary solution of (25). 

Given a surface distribution of sources, it can be shown that the com- 
ponents of the gradient and therefore of the hyperbolic gradient of the 
potential remain finite and continuous on either side of the surface S. 
Also, ®, and therefore its tangential derivatives, are continuous across the 
surface. 

In order to find the discontinuity of the normal derivative across S, we 
apply (41) to a small cylinder whose bases are parallel to the surface on 
either side of it, and whose height is again small compared with its lateral 
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dimensions (Fig. 2). Letting first the height of the cylinder tend to 0 we 
find that 

[ gradh®.dS— { gradh®.dS = 2z [ od, 

Ss, £.. s’ 
where S’ denotes the portion of S inside the cylinder, and S’, and §'. 
denote the two bases of the cylinder respectively, S',. being the base whose 
outside normal coincides with the outside normal of 8. Letting S’ tend 
to 0 round any given point on S, we obtain, denoting by A, yp, v the 


h, 





a P) 
NX 
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Fia. 2. 


direction cosines of S at the point in question, and indicating by + the 
derivatives of ® on the two sides respectively, 
— Ne nl) oP) — ano, 
Ox, Om. Cy,  ey- Gz, O2_ 
In particular, if the distribution is in the (x, y)-plane, we have A = p = 0, 
v = 1, so that 





cD af ‘ 
— ——_ = 270, (45) 
CZ 4. CZ 


a result which is of considerable importance for the calculation of the 
wave drag of an aerofoil moving at supersonic speed at zero incidence 
(refs. 4, 5, 6). A similar relation for the discontinuity of the potential 
across a doublet distribution in the (x, y)-plane, and which can be derived 
from (42), is fundamental in the supersonic theory of flat aerofoils at 
incidence. 

These relations have hitherto been inferred by analogy with incom- 
pressible theory and then proved ad hoc in the particular cases required. 
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The need for a more systematic development was pointed out in the 
introduction to ref. 5. 

Dividing (44) by A = ,/{(A8)?+p?+v*} we obtain the result that the 
discontinuity of 0®/én’ in a direction n’ whose components (direction 
cosines) are (- ¥ x) is none And since the tangential derivatives 

AAA A 
of ® are continuous across the surface, it follows that the discontinuity 
of @D/2n’ must be the discontinuity of 0O/én, where n = (A, p,v) is normal 
to S, multiplied by the cosine between n and n’. Hence 





& jet = 2mo0(—A*B?+-p?+v?)/A 

en on 

co o@® , 1—A*(1+-f?) 
) Senet anon - = 2 ———— eae e 46 
7” on, on Jfl—r2(0— p)} aie 


Equation (46) amends the statement in ref. 5 that the discontinuity 
of the normal derivatives is always 270. The particular case in which 
this statement was applied, however, viz. (45), remains correct. 

The chief use to which Hadamard puts his concept of the finite part is 
related to the above applications but is the outcome of a rather different 
approach. Hadamard’s purpose is the solution of Cauchy’s initial value 
problem for a very general class of hyperbolic partial differential equations 
including (20) as a special case. We are going to develop Hadamard’s 
result in respect of equation (20), i.e. we are going to find an expression 
for the value of a solution ® of (20) at a point P = (29, ¥9,2») inside a 
closed surface S, when the values of ® and of @@/én’ are known on S, 
where, for every point of S, the direction n’ is defined as above. An 
equivalent formula had been derived previously by Volterra (8). 

Let ‘¥'$ be the velocity potential of a counter-source of unit strength 
located at P; then the function ®Y} is admissible inside S, excluding 
only P, provided ® is regular inside S—and this can in fact be verified 
a posteriori. 

Let us surround P by a small surface S’ and apply equation (32) to the 
volume V bounded by S+S’. Since Ah® = Ah} = 0, we obtain 

° 
(® gradh ¥$—YV$ gradh®).dS = 0, 
S+S’ 
or, taking the inward normal as the direction of a surface element of S, 
and the outward normal as the direction of a surface element of S’, 


* * 
| ®gradh ¥$—¥V$ gradh®).dS = { (®gradh ¥$—¥$ gradh®).dS. 


(47) 
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It can be shown that as S’ contracts to the point P, the left-hand side 
of (47) tends to 27®(xp, yo, 2%). Thus 
* 
2rD(X, Yo: %) = | (@gradh ¥$ —V3 gradh®).dS. (48) 
§ 
Denoting by A, p, v the direction cosines of the normal to dS, we have, 
for any scalar function ®, 
gradh®.dS = (- ye use Sas. 


a 


For 8 = 1, this becomes gradh®.dS = -*. dS, where n* = (A, —p, —»). 
n 
Hence, for 8B = 1, 


* 
at) - ds. 


2nD(2x9, Yo: %) = | (we r (49) 
FS 


n* — @n* 
This is Hadamard’s formula (58) (ref. 1, p. 207) for the special case 
f = 0. The direction n* is called by Hadamard the transversal direction 
to dS. Its geometrical interpretation, due to Coulon (compare ref. 1), is 
that it is conjugate to the tangent plane to dS with respect to the cone 
(w—a,)?—[(y—y,)?+(z—z,)?] = 0, whose vertex is located at dS. 


5. Vorticity distributions in steady supersonic flow 


We now direct our attention to the study of rotational motion. 

Given a field vector q = (u,v, w) we denote by &, n, ¢ the components of 
curlq, € = al 9 = ee, The differential equation 

Cy OZ cz Ox ox oy 

of the system of vortex lines is dx/€ = dy/n = dz/{ as usual, and the 
strength of a vortex tube is defined as the product of the cross-section ¢ 
into the resultant vorticity w = (£*+ 7?+ C?)!, and is the same at all points 
of a vortex. All these results and definitions are in fact quite independent 
of whether the fluid is compressible or incompressible, except that in the 
case of supersonic flow it may be necessary to consider the finite parts 


of integrals of the type { q.dl and f curlq.dS in cases where the ordinary 
C 8 
integrals do not exist. 

Applying the vector operator Vh in cross-multiplication to a vector 
q = (u,v,w), we obtain a vector which will be called curlh q (hyperbolic 
curl of q). 

Explicitly— 
wthe ("- év ou ow dv “| (50) 


ee a | 
ey az’ dz +P a “ oy 
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Direct calculation shows that 


divh curlhg = 0 (51) 
and curlcurlhg = gradh div q—div gradhq. (52) 


A field vector q will be called irrotational or lamellar, as usual, if 
curlq = 0, and it will be called hyperbolic solenoidal if divhq = 0. 

We are going to show that a vector q defined in a region R and admissible 
init can be represented as the sum of three vectors, one irrotational, one 
hyperbolic solenoidal, and one both irrotational and hyperbolic solenoidal. 

More precisely, we are going to represent q as q = q,+q,+q,3, where 


divhq, = divhq, curlq,=0 in R, (53) 
curlq, = curlq, divhq,=0 in R, (54) 
and divhq, = 0, curlq, = 0 in R. (55) 


Assuming that vectors as described in (53) and (54) have been found, 
we put q,; = —q,—q,. Then divhq, = divhq—divhq,—divhq, = 0, 
and curlq, = curlq—curlq,—curlq, = 0, so that q, defined in that way 
satisfies (55). 


, Ds ; F , 
Putting o = —divhq in R, we determine a scalar function ® by 


“a7 
* 
=| oDEdV,so that Ah® = —2z7e according to (42), ic. Ah® = —divhq, 
R 
and so q, = —grad® satisfies (53). Thus 
* 
q, = — grad [ divng.of dV. (56) 
pi R 


To find q,, we shall assume that q, is given as the hyperbolic curl of 
avector ¥ = (‘f;, ‘V4, ‘¥3), dg = curlh¥, and so by (52) 





curlq, = curlcurlhY = gradh div ¥—div gradh®. (57) 
We now restrict Y by the condition that div'¥ — 0. Then we must have 
div gradh ¥ = —curlq, = —curlq, by (54) and (57), ice. 
Ah¥Y, = —£, Ah¥Y=-—n, Ah¥,=—C in R, (58) 
and this according to (42) is solved by 
S ba * 
¥%,=1([gsav, v= [ rosdV, Yat | cox dV, 
2m, + 27 J 2a 
R R R 
ca 


= y — = | curlq®} dV. 


R 
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Then q, = curlh¥ satisfies (54), provided we can show that in fact 
div'¥ = 0, as assumed. And this can be shown exactly as in the classical 
counterpart (ref. 7, para. 148), provided the integrand vanishes at infinity 
or, in particular, if it vanishes outside a finite region. This again will 
certainly be the case if curl q vanishes for sufficiently small x, since ©, 
vanishes for sufficiently large x. 

It is clear from the above construction that q, determines the flow due 
to the source distribution in R, while q, represents the flow due to the 
vorticity distribution. Given a three-dimensional vorticity distribution 
we then have in detail 


* 


dx, dy, dz 
W(x, y,z) = x | | (é, ; ¢) V[(w—a)? ri — 0 + (2— Z)*}] (60) 


(the expression on the right-hand side is actually an ordinary improper 
integral), where R’ is the sub-domain of R which satisfies 





FR’ 


(w—2ao)? > B*[(y—Yo)? + (z—2)?] and ay <2. 


We then obtain for the components wu, v, w, of qg = curlh'¥, 














* 
p2 | : dx, dy, dz 
w = — | [(2—2%)n—(y—Yo)£] By —y,)2 zy 
al . [(x —B{(y Yo)" + (z—Zo)* }] 
a: dx, dy, dz 
= 5. | [(x—a)b—(2 lia) — BX a Pte 2o)*}}* (61) 
ge 
2 dx, dy, dz 
v= — * ftw—we- (x e— to) ae \ — Bf ( aie z—Z2,)*}]} 
" 





This may be written 
* 


, 
qd. = — [ (rA curl) , (62) 
& 
where r = (w—2, y—Yo, 2—2), and s = [(a—ay)?—B*f(y—yp)? + (z—20)}}* 
The corresponding formula for incompressible flow is 
* 
q, = [ (r A curlq’) il (63) 
R 


The discrepancy in sign is only apparent, since as, for instance, formula 
(8) shows, the sign of a finite part does not follow the sign of the integrand, 
as for ordinary integrals. 

We may now calculate the field of flow due to an isolated re-entrant 
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jine vortex C. Replacing the volume element dz, dy, dzy in (60) by ody, 
where dl, is the element of length of C, and o, its infinitesimal cross- 


ection, and writing w = (€?+ 7?+?)!, we have 
_ day — yy {=o deg 
; dl, dl, ll, 


and so, since wo, is a constant, K, and since dl, = (dx, dy, dz), 


60) becomes 
* 


K f dl, : _ 
2a J yi (x - Xq)®—B*{(y—Yo)? + (z—Z9)}]’ 


( 


Wax. Y,Z) 


a 


(64) 


where C’ consists of the segments of C which satisfy 
B*|(y—yo)?+(z—Z»)?] and 2x > ay. 


Ifin particular C consists of straight segments which are parallel either 
to the 2-axis or to the y-axis, then (64) can be integrated, since 





| dry ) 
) V[(w—a9)?—Bf(y—Yyo)® + (z —2)*} ] 
x—x 
cosh-! —— =o = + const. 
Py(Y—Yo)" + (@— Zo)" 
and r+ (65) 
1 dy 
. V[(a Xo)” Bt(y Yo) \< Zo)? | 
Diss y—2 
—-—sin-! 3 es Yo) 5, + const. 
ad V{(®—29)>— B?(z—2)?5 J 
Now let C be a horseshoe vortex of strength K, consisting of the straight 
segments (7, < %y <D, Yo Ya» %q = 0), (Zo = 2%, —Hs S Ho < Has % = O) 
and (x, < %y <<, Yo = Y;, % = 0), where x, and y, are given constants 
Fig. 2). Using (65), we find that ‘Y, = 0 always, and ¥, = ¥, = 0 for 
t <2, while for x > 2,, 
: K x—xX x—2X 
7, — | cosh-! —— L > — cosh aT, 1 —|, (66) 
“7 Py Uy Yy)° +2; By, UY+Y1)° +275 


where the cosh-! are to be replaced by 0 when their arguments are smaller 
than 1, respectively, and 





' : 3(y—yYy,) , (y+2 ies 
Y sin-! A —sin-! - _Ply+yn) _ | (67) 
2 3.2 T= » \2__ Q2,2) Tins ae 
27D (a x) Bz?) i (a x1) B22? 
Where the sin-! are replaced by +42 or —47 when their arguments are 
greater than 1, or smaller than —1, respectively. 


5092.4 
Ff 
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We now obtain q, = (u,v,w) by taking the hyperbolic curl of 
Y = (4, 4,0), so that w = v = w = 0 for x, > 2 while for x > z,, 
— «lr NE ee Z 
2a |[(x—2,)*— 2? ]{(x—ax,)?— BL (y—y,)?+2"]}! 





4 


to! 


~ [(e—ay)P— Re] x—2,)? i [( (y+y,)? 


t 





7 (x—2z,)z 
~ [yt+y)?+22]{(@—2,)?— Gah 
=| _(®@—a)(y—y ){(w—% By) P+ 227]} 
[(e—a,)?— 2" [(y—y,)?# +2 {(e—a,)?-BP[(y—yy 2 +2} 
W—Bytysyr+ 220] 
*| y,)?-+2]}! i 











(x—2, My tw tl: 


nom 
~ [(@@—a,)2— Belo ty Pra 


\(x—a)?#— B? (ys 





(68) 


where, for given x, y, z, the imaginary terms are omitted. Except for the | 


notation, equations (68) agree with the field of flow round a horseshoe 
vortex calculated by Schlichting by an entirely different method (ref. 2). 

Some care is required when attempting to represent a volume or surface 
distribution of vortices as a combination of line vortices. Thus, according 
to (68), the ee u, v, w, all vanish when (x,y,z) is outside both 
the cones (w—2,)?—f*|(y+y,)?+2?| emanating from the tips. But it can 
be shown that this is no longer the case if the vorticity in the spanwise 
segment is distributed over a finite width Az). However, even then the 
failure (which is due to the discontinuity of ‘’ for line vortices) can occur 
only at points belonging to the envelope of the cones of type 


. - \2 2 -. i 
(x—2)*—B l(y— Yo)” +(z—Zp) ] = 0 
emanating from the vortex lines which are supposed to generate the surface 
or volume. 


6. Review of principal points 
The present paper contains various developments of the analogy 
between linearized supersonic flow and classical incompressible hydro- 
dynamics, with a view to applications in supersonic aerofoil theory. 
The equation satisfied by the velocity potential in linearized supersonic 
flow is 
pe a @, (69) 


on Of ' 
where 6? = M?—1 = U*/a*—1, M = U/a being the Mach number, 4 the 














yeloc 
direc 
direc’ 

Elk 
defin 


and 


whel 
pote 
pote 
in a 
—U. 

Tl 
by ¢ 
of J 
doul 
who 
equi 
that 

T 
in § 
as 0 
soni 
sucl 
and 
in § 
(on 
The 
dist 


late 


for 
Fir 
floy 
cal 
wil 
the 


url 0! 


(68 


or the 
seshoe 
ef. 2 
urface 
rding 
> both 
it can 
unwise 
on the 


occul 


urface 


valogy 
1ydro- 


rsonit 


(69 


a the 





f 








LINEARIZED THEORY OF STEADY SUPERSONIC FLOW 431 





velocity of sound, and U the magnitude of the free stream velocity. The 
direction of the free stream is supposed to be parallel to the positive 
direction of the x-axis. 
Elementary solutions of equation (69) are the functions ®,(z, y,z) 
defined by 
D>(x, y, 2) SS 
| (x—2%q)?—B*{(y—Yo)? + (2 —Zo)*}] 70 
for (x—xy)* > Bl (y—Yo)?+(z—2)*], x > Xo (70) 


and @,(x, y, 2) 0 elsewhere. 


where P = (Xp, Yo,%) and o are arbitrary. ®, is said to be the velocity 
potential of a source of strength o located at P. The actual velocity 
potential of a physical (acoustic) source travelling with a velocity —U, 
ina frame of reference travelling with the source is obtained by adding 
—Ux to Op. 

The velocity potential of a supersonic doublet is, by definition, obtained 
by differentiating ®, with respect to length, relative to the coordinates 
of P in any given direction. As in incompressible flow, a supersonic 
doublet can be interpreted as a combination of two supersonic sources 
whose distance 6 from each other tends to 0, and whose strengths are of 
equal magnitude and of opposite sign, +o’, and related to 6 in such a way 
that the product o’6 remains finite as 5 tends to 0. 

The importance of these and other particular solutions of (69) mentioned 
in §§ 4 and 5 above depends not so much on their direct physical meaning 
as on the fact that the flow round aerofoils and other bodies under super- 
sonic conditions can be built up from these solutions. With a view to 
such applications, general line, surface, and volume distributions of sources 
and doublets are considered in §4, and similar distributions of vorticity 
in$5. In §4 results are derived corresponding to the theorems of Gauss 
on total normal intensity) and of Poisson in classical potential theory. 
The discontinuity of the normal component of velocity across a surface 
distribution of sources is calculated as an application. 

The velocity field due to an arbitrary distribution of vorticity is calcu- 
lated in §5. In particular, formulae are derived for the field of flow due 
to an isolated re-entrant vortex travelling at a steady speed. These 
formulae constitute a supersonic counterpart to the law of Biot-Savart. 
Finally, it is shown that for the case of a horseshoe vortex the field of 
flow obtained by the method of the present paper agrees with that 
calculated by Schlichting in an entirely different manner. In contrast 
with the fictitious character of sources and sinks in supersonic aerofoil 
theory, the occurrence of trailing vortices in the wake of an aerofoil is 
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a matter of physical reality. On the other hand, the vorticity induced 
by a shock wave is not taken into account by linearized theory. 

Two mathematical ideas are required to carry out the above analysis, 
One is the purely formal introduction of the vector operator Vh £ (‘hyper- 
bolic nabla of index 8’) defined by (—# A meh which is used in 

ou Cy &z 
conjunction with the ordinary vector operator V = (; a = The 
Ox’ Cy Oz 
hyperbolic nabla, like the ordinary nabla, can be applied either to a scalar, 
or, in scalar and in vector multiplication, to a vector. 

The second concept required for the analysis of the present paper is 
that of the finite part of an infinite integral. It was introduced by 
Hadamard in connexion with the solution of Cauchy’s problem in partial 
differential equations (ref. 1). This concept enables us to cope with the 
infinite integrals that occur in the analysis of supersonic flow, e.g. in 
calculating the flow across a surface surrounding an isolated supersonic 
source. For the purposes of the present paper, the concept had to be 
developed in rather greater generality than considered by Hadamard. 
Also, in addition to Green’s theorem for finite parts, a particular case of 
which is proved and applied by Hadamard, we require a proof of the 
extension of Stokes’ theorem to include finite parts. 

The analysis of the present paper has been applied to the calculation 
of the downwash in the wake of a Delta wing. The results of this investiga- 
tion are given elsewhere (ref. 9). 
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THE APPLICATION OF THE NATIONAL ACCOUNTING 
MACHINE TO THE SOLUTION OF FIRST-ORDER 
DIFFERENTIAL EQUATIONS 
By A. E. CARTER and D. H. SADLER 
H.M. Nautical Almanac Office, c/o Royal Observatory, Greenwich, S.E. 10) 


[Received 2 December 1947] 


SUMMARY 
Milne’s formula for approximate quadrature is used as the basis of a method for 
the solution of first-order differential equations on the National machine. The 
method, which is illustrated by a numerical example, enables the machine to form 
the required dependent variable without the necessity for conversion from a sum 


to an integral. 


1, Introduction 

Tur National Accounting Machine is now in fairly general use for scientific 
computing; it has been described in detail by Dr. L. J. Comrie (1), to 
whom belongs the credit for realizing the possibilities of the machine and 
for exploiting them. For those to whom the machine is unfamiliar, a very 
brief description of its principal features is given in the next section. 

The machine has for long been applied to assist in the solution of 
differential equations of both first and second orders, but its main use has 
generally been restricted to saving the operator the labour of writing and of 
differencing and integration. Since the machine has a limited storage 
capacity and there is no means of direct multiplication, there are obvious 
difficulties in using the machine to convert derivatives to the differences 
which are necessary for building up the solution. A method for the solution 
of linear second-order differential equations, which overcomes this diffi- 
culty, has already been described (2); this method relies, however, on a 
mathematical transformation of the original variables. 

The solution of first-order equations suffers from the disadvantage that 
first derivatives become available at the tabular interval, while the first 
differences are required at the intermediate half-way point. The method 
given below, in its simplest form, relieves the operator of all stages other 
than the calculation of the first derivative, which cannot be avoided. 

2. The National machine 

The National is a printing adding machine with six adding mechanisms, 
or registers, each with 12-figure capacity; two of these registers (Nos. 1 
and 3) can subtract directly as well as add, but subtraction in the other 
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four must be performed by the addition of complements. A number set 
on the 12-bank keyboard can be added (including subtraction in registers 
1 and 3) to any combination of the six registers; the contents of any 
selected register can be transferred, with (a total or T operation) or 
without (a sub-total or ST operation) clearing of the register selected; in 
each case the number added or transferred is printed. The destination 
of such a number is controlled by the ‘stop’ on which the machine is 
operated, and thus, in any particular set-up, by the position of the moving 
carriage; register selection can be made by the operator independently 
of the carriage position, but is in fact almost invariably governed by 
position. 

The moving carriage carries a removable ‘form-bar’ on which stops 
can be easily placed in any desired position and order; it normally ‘tabu- 
lates’, or moves from one stop to the next, after each operation, automati- 
cally feeding the paper and returning to the first position, indicated by 
an M (margin) stop, after reaching the end of a line, indicated by an R 
(return carriage) stop. 

It will suffice here to mention further four points of detail; stops can 
be made available to add in any combination of registers and, in particular, 
‘reversible’ stops are available for registers 1 and 3 enabling them to be 
changed quickly from + to — and vice versa; the machine can be made 
to interpret the instruction to subtract the contents of a register from 
itself as an instruction to clear the register concerned; printing can be 
suppressed in any position by using an NP (non-print) stop; an R stop can, 
when necessary, be rendered ineffective to allow of the carriage moving 
to positions beyond it. 


3. Notation and basic formulae 
The differential equation will be written in the form: 
di 
oo f(x, y). (1) 
dx 

The interval of tabulation is denoted by h, but, as is usual, consecutive 
values of the independent variable will be indicated by integers and 
corresponding values of y and f by suffixes. Central difference notation 
is used. 

The fundamental requirement is for an equation connecting two, not 
necessarily consecutive, values of y with the derivatives, and their 
differences, at intermediate points. If the coefficients of the derivatives 
and differences can be reduced to small integers it may be possible to 
perform the simple multiplications and additions on the machine. 
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APPLICATION OF THE NATIONAL ACCOUNTING MACHINE 


The following formulae are available: 


Yn—Yn-2 2h( f+ 48?— 7484+...) na 

Lh(f,,-o+4f,.+f,,)—hd'+... (Simpson’s rule) ; (2) 
Yn—Yn—4 4h( f+ 387+ 764— 57,6°+...),-» 

th(2f,,-3—fn—-2 + 2f,-1) +){h84—... (Milne’s formula); (3) 
Yn—JYn-6 6h( f+ 35? T 50" T ai50° — 36590 ° +--+) n—3- (4) 


The formulae with an odd number of intervals are generally unsuitable. 
Of those given (3) is by far the most tractable. Putting 


g(x, y) = $hf(x, y), (5) 
Yn—Yn 4 (3g+ 287+ 54 —47,0° +...) no, (6) 


where the differences now refer to g. The formula is used in this form 
rather than in the Lagrangian form for a number of reasons. Firstly, the 
National machine lends itself more readily to a formula containing 
differences instead of a series of ordinates, secondly, differences give some 
form of check on the integration, and thirdly, the fourth differences are 
required in order to allow for their effect. 


4, National set-up 

The principle of the set-up is to difference g, to the fourth difference, 
at the same time building up (3g+26?),_, in an empty register; the 
addition of y,_, will then give y,..,, which is used to calculate g,,, and 
thus to reeommence the cycle. Provision is made for applying the fourth 
difference correction without having to reset the machine. 

The set-up is given in detail in Fig. 1; much of the complication arises 
from the necessity to cater for varying signs. In the following description 
such details will be ignored; these can be followed from the set-up by those 
possessing a National machine. 

In addition to the National, a calculating machine is required for com- 
puting g from a knowledge of x and y, and a second machine or table for 
Inding 784, 

The normal cycle consists of: 

(i) Operate the machine in positions 1, 2 (printing the argument and 
corresponding approximate function), 3, 4 and 5. 

(ii) Caleulate g, from the value of y, printed in position 2; in doing 
this arrange, if possible, that y is the last multiplier or that there 
is some simple method of correcting g differentially for a small 
change in y. 

(iii 


Set this approximate value of g, in position 6 and operate the 
machine in position 7 to print the corresponding 64 _.. 
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APPLICATION OF THE NATIONAL ACCOUNTING MACHINE = § 437 


(iv) Using the table, or second machine, calculate 764_ 


2 and use the 
difference between this and the extrapolated value (appearing in 
position 14 on the previous line) to correct the approximate value 
4 


of y,,; change g,, accordingly, and 54_, by the same amount; if this 


change is large, it may change 464_., y,. g, and 54_ 


2 itself; but 
the process is rapidly convergent and a balance is soon achieved. 
In practice, the interval h will be chosen so that the extrapolation 
can be done with little error. The final correction to y,, must be 
recorded, by hand or otherwise, either directly or in the form of 
an alteration of y,, and the final 64_, entered in position 8. 


< 


Operate the machine in positions 9, 10, 11, 12; these, in conjunc- 
tion with the operations in positions 3, 4 and 5, result in putting 
3g-+ 28 into register 3, as well as completing the differencing cycle 
(vi) Set y,,, from the corrected value in position 2 on line n—3. 


(vii) Extrapolate 54_, and set ,464_, in position 14. 

At convenient intervals, say, every 10 or 20 lines, the accuracy of the 
various settings should be checked and 84 should be differenced for 
checking purposes; a check on the accuracy of the machine is provided by 
comparing g,_, printed in position 5 with the value printed in position 6 
on the line above, corrected for the difference between the approximate 
and final values of 54_,. The set-up suffers from the disadvantage that y,, 
is checked only through the differencing of g,,; therefore the contributions 
to y, set into the machine in positions 13 and 14 should be checked, since 
small errors might otherwise be undetected. 

It might be thought that this cycle, requiring operation of the National 
machine, an ordinary calculating machine, reference to tables and some 
mental arithmetic, would impose a considerable strain on the computer. 
Actually, all the operations other than the extrapolation of 54 and the sub- 
sequent correction of y soon become almost automatic and the computer 
can concentrate upon the general flow of the solution, thus achieving a much 
higher accuracy in extrapolation than would be the case if his attention 
had to be diverted by writing or by mental additions and subtractions. 


5. Accuracy of the method 

The method can be made as accurate as desired by retaining sufficient 
figures and, if necessary, decreasing the interval. Although no provision 
is made in the set-up, it is not difficult to allow for the sixth difference 
correction by extrapolating it and setting after position 14; but if 8° is 
large enough for ,3;5° to be sensible, extrapolation of 54 will have become 
difficult and either a change of interval or a reduction in the number of 
figures retained would be called for. 
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The building-up error, which cannot be avoided in numerical integration, 
is here quite different from that normally occurring. There are really four 
separate series, of which a typical one contains the sum of errors of the 


form: 
. 9 Lé¢ 5 
Z€n—gt€n-2t 2€n-1 + Nn-2) (7) 


where the e’s and 7’s are elementary rounding-off errors. This compares 
with the straight sum of four rounding-off errors, which is the optimum 
that can be achieved for four steps. The building-up error is therefore not 
likely to be serious and the capacity of the National will allow ample 
guarding figures to be retained. The errors of the different series are not 
independent so that they cannot be determined by examining the 
differences; the second-order effect of these errors on the solution is beyond 
the scope of this paper. 


6. Example 
To illustrate the principles of the method it will be applied to the 
equation es 
w= aty 8) 
dx 


with initial condition y = 1 at x = 0. This example has been deliberately 
chosen because of its simplicity and the constancy of sign of the differ- 
ences. The interval / will be taken as 0:1, so that 
g = 0-133...~+-0°133...y. (9) 
Nine decimals will be retained. 
By expansion in ascending series, by successive approximation, or by 
other methods (e.g. repeated use of Taylor’s theorem) we find: 











Line eS 4 y g and differences (unit gth decimal) 
-2 —o'2 +0°83746 1506 + 8499 4867 
| +2296 1778 
—I —-o1 | 0:90967 4836 10795 6645 +241 4910 
2537 6688 +25 3980 
° oo | I*00000 0000 13333 3333 266 8890 +2 6710 
2804 5578 28 o6g0 
I +0or1 111034 1836 16137 8911 294 9580 
3099 5158 
2 +0°2 1°24280 5516 19237 4069 





As a check it is noted that 
(3g+ 26+ 54), = +0-40534 4011, 
whereas Yo—Y-2 = +0°40534 4010. 
These initial values are fed into the National according to the rules 


given on the set-up and the computer is faced with the figures appearing 
above the line in Fig. 2. Note that the value at 0-2 differs from the 
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computed value. The next stage is to extrapolate {5 and the computer, 
very much in the dark, chooses the convenient round number, 7000, 
which is set in position 14. 

Stage (i) is then completed and the computer turns to the calculating 
machine to form g, from the recently printed y,. Here 0-133... is put into 
the product register and 0-133..., on the keyboard, is multiplied by the 
printed value of y, to give an approximation to g, which is set in position 6, 
The next step is to print the approximate fourth difference 5{, namely, 
2 9534; this shows that 7000 was a very poor extrapolation for #84, the 
true value being 6891. This change enters directly into y, and therefore 
into the multiplier on the calculating machine, which is now reduced by 
109 = 7000—6891; this changes g, by 15, 5? by the same amount, and 
thus 6} and y, by 3; g, is, however, unchanged and so is 3457. The whole 
process is shown in tabular form below: 








> >R s4 _7 94 

M.R. = y3 P.R. = g, oj 3991 
Extrapolated value 7000 
First approx. 1°3997I 7729 22662 9031 2 9534 6891 
Second approx. 7620 go16 2 9519 6888 
Final values 7617 go16 





The value of y, is now amended by hand (or the correction to 744 and 
Y3 printed by the machine) and the final value of 54 = 2 9519 is set in 
position 8; there is no need to correct g, since the correct value will be 
built up from 6}. Routine operation then brings the machine to position 13 
in which the corrected value of y, (here 1-00000 0000) is carefully set. In 
the next position 453 must be extrapolated from the two printed fourth 
differences; a value of 7600 is chosen and the cycle recurs. 

It will be noticed that, after the first, the extrapolated values are all 
within 14 of the final values, leading to a negligible second-order effect. 

After six more cycles we obtain: 


at x= 1:0, Yy) = 3-43656 3668 as compared with the 
true value = 3°43656 3657. 


Now it will have been clear that sixth differences, although significant, 
have been ignored; since their effect is very small, the correction need not 
be applied to, say, y, until this value is required on line 6, by which time 
5$ is accurately known. The application of the two relevant corrections 


2.58 = 296 __ 
—s3ipo4 = —3 and —,;7;08 = —4 


reduces the discrepancy in the check value. But even when allowance 1s 
made for the error of —1 in y, there is still a difference of 3, which is 
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rather large; examination shows that this is due to an unfortunate series 
of rounding-offs in g. 

Detailed modifications are clearly possible; in particular the final correc- 
tion tu y can be set on and printed by the machine instead of altering the 
end figures of y by hand. 
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ON THE HODOGRAPH TRANSFORMATION FOR 
HIGH-SPEED FLOW 
II. A FLOW WITH CIRCULATION 
By M. J. LIGHTHILL 
(Department of Mathematics, The University, Manchester) 


[Received 13 December 1947] 


SUMMARY 
A steady two-dimensional isentropic compressible fluid flow about a contour with 
circulation, reducing to the potential flow with circulation about a circle as the 
Mach number at infinity tends to zero, is specified in the hodograph plane. 


1. Introduction 

In this paper the method of Part I (1) is applied to obtain the high-speed 
flow corresponding to the low-speed flow with circulation about a circular 
cylinder. The reasons for presenting this alternative treatment of the 
problem to those of refs. (2) and (3) are again those attaching to any such 
multiple treatment: it widens the basis of understanding and future 
development. The work of this part is necessarily more complicated 
than that of Part I, however; nor has it the generality and directness in 
argument of ref. (2). Its mathematical interest, in the author’s belief, 
outweighs these objections. 


2. The low-speed flow 
If the circulation is 47sina and the velocity at infinity unity, the 
complex potential and complex velocity for the flow of ‘perfect’ fluid 
past the unit circle are 
l — , dw 1  2isina 
w = 2+-+4+ 2isinalogz, (2 6 (1) 
z ; dz 


ae ~ 


_ ~ 


Inverting the latter equation as z = (1—{)-1(—isina+(cos*a—{)#), we 
obtain for w the expression 
1 _ ce 
w = —~(—isina+(cos*a—Z)*) +7 sin a+(cos*a—l)!— 


| 
—2isin «log(isin a+(cos*a—{)?) = P(f)+Q(g), (2) 
where P({) = isin a(1— ;— lost), (3) 


and Q(¢) can be put into the form 


Q(¢) = (1+ = _)tt—t—sinta)!+2 in vsin-"((1—{)-!sina). (4) 
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By expansion of @ in descending powers of (1—{), which are in turn 
expanded in ascending powers of ¢, and by inversion of this double sum, 
we find that the branch of Q(Z) equal to 2(cosa+asina) when = 0 has 
for its Taylor series near £ = 0 the expression 


= be (n 1)(n—#)! 


an n! (—3)! 


F(n—4, —4; 4; sin*a). (5) 
The hypergeometric function appearing here is an integral function of n: 
it has the integral representation (for all 7) 
sin?a 
(cos «)!-?”— (n—4)sin « a-*(1—a)-!-" dz, (6) 
i) 
which shows that it is O( n(cos*«)-"|) as |n| +00 for R(n) > 0. 
Hence the series (5) converges for |f| < cos*a. When R(n) < 0 expres- 
sion (6) becomes 
; 
—(n—4)sina | x-#(1—ax)-!-”" dx+O(|n(cos?ax)-"|) 


0 


— se)! (—t-d)! . 
. (—3):- ( . 3) sin a+ O(|n(cos*a)-"|) = O(|n|). (7) 


( nN): 


The series for Q(f) can be continued in the manner adopted in §2 of 
t > 
Part I. In fact the series is 


—8)' Fw—4, —454;sin*a)(—L)" dv (8) 


within its circle of convergence |f| < cos*a, being minus the sum of the 
residues of the integrand at its poles y = n > 0 to the right of the contour; 
while for || > 1 the integral, by (7), is equal to the sum of the residues 
of the integrand at its poles to the left of the contour, i.e. to 


3)! 


vi (—2—b(2— 


— —)>, 
n=0 ao 


F(—n, —4; 4; sin’a)(—1)"(—f)#-. (9) 


y 


Here, as in Part I, if argf > 0, arg(—) must be taken as arg{—z in 


order to secure the convergence of the integral (8), so that (—{)!-" is 


—i(—1)"f-"; while if arg < 0 then arg(—Z) is arg{+-7 and (—{£)!- is 
(—1)"-". Thus, for |f| > 1, if + signifies the sign of arg on the path 


of continuation. 


, _. <> (—n—4)(n—3)! eis 
Q(¢) =| > 1—3n—3) F(—n, —4; 4; sin*a)f?-". (10) 


n!(—4)! 
n=0 - 
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This can be continued into the region cos*a < |{| < 1 by use of (7), whence 


Qo) = Fr > = — : Fe gin! sin a+ O(n cos?" »)| 2 n 


n!(—+4)! n—+4)! 
oe (—3) v—$)! 
x 
1 
fies tin. 
= Fi > =" sina f?-” (11) 
L_yp 
Dy) 
n=0 “ 


plus a function analytic for |f| > cos*a; (11) in turn differs by another 
function analytic for |{| > cos*« from 

isin af (1—f)-!+log(1—Z)], (12) 
with which equation (3) should be compared. It is deduced that 
P(f)+Q(¢), continued analytically via the region arg f < 0, has no singu- 
larity at ¢ = 1; but that if continued via the region arg > 0, it has a 
singularity there, that of —2isino{(1—{)-1+log(1—Z)]. Also, by (10), 
Q(¢) changes sign in going round a cut joining cos*x and 1. The Riemann 
surface corresponding to that in §2 of Part I can now be depicted as 
below, Fig. 1. The values of w shown are those obtained by continuation 


arg 5-0 





Fic. 1. 


along paths not crossing a cut from infinity to the singularity S: w is 
one-valued on the surface so obtained. AB is the cut joining the two 
sheets. A positive encirclement of S adds 47sina to w. The only singu- 
larities of w are S and B, corresponding in the physical plane to the point 
at infinity and a certain finite point C where z = —icoseca. The cut AB 
corresponds to a circle of diameter OC (O the origin) which passes through 
the two stagnation points (corresponding to A): on this circle 6 = 0. 


3. The high-speed flow with f,,(7,) = e-" 
To obtain from this a high-speed flow, we write, as in Part I, § 3, corre- 
sponding to the above P, Q, 


P= —isina x 1—n-)us,, (r)e~ "81 +#9), (13) 


—1)(n—3)! ‘ ;, ; 
Q= > CoO D! raat: dssintaypa (rien, (18) 
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and consider 3(P+Q) as a commencing branch of y%. The former series 
converges for + <7,, the latter (by equation (6) and the asymptotic 
formula for &,,(7)) for e* < e cos*@: since s is an increasing function of 7, 


this can be written + < 7, for some 75. 
When 7 < 7, we have as in (8) 
] r Gs 1 )( v Ll)! (yp i)! a we . 
Q= — =!* F(v—4, —4; 4; sin2a)pb, (r)e"F i711 dy, 
271 (—4) - 
(15) 
with the + sign according as 6S 0. For rt > 7, this is equal to the sum 


of the residues to the left of the contour, i.e. 


_ ~ (—n—4)(n sy! , _ 
q=7F:1)5 st _2"* F(—n, —}4; 4; sin2a)ypy_,, (rer —Der t+ 
ra nN: | -})! . 
nm 0 


« (—n—1)(n—1)! (—n—3)! 5, -— 
- > 2) F(—n—}, —}4; 4; sin*a) x 


< (—n€,, p,(7))(—1)"er™+, (16) 

The latter series, which we call R, is convergent for all 7, 6. Therefore 

the continuation of Q round a cut from 7 = 7,, 0= 0 tor=7),0=0 
is 2R—Q; and continuation a second time round gives again Q. 

Also, continuing the series for -(Q—R) in (16) into the region 


T 7 71, we obtain 
(Q— R) 
1 ;\? 1 ' ' 
, < n 5 )(n yt Sin: . ‘ a 
iS 5 kgf ee sin «+ O(n cos?"x) | by_ (ren —tXei +18) 
1)! 1\! rs—n\ 
pa nr: | y): | (m—s): 
? 0 - . 
l 
x" = mw ., 1 . ~ 
> a sin a yy _ (rel X81 +00) (17) 
l 2 
— 5 Mt 
0 = 
plus a function analytic for 7 >7,. Expression (17), for 7 > 7,, is equal 
since the residue of z/coszv at v 1» is (—1)"-") to 
sin a : 7 v—|] (tin—s—i0 - ' 
— ob, (ren Fir—si1-1) dy +X (18) 
271 F COS 7TV Vv 
vccording as @ < 0, where 
ba sin a s a(n+1)C, us,,(7)¢ n(81 +18) (19) 


0 
“9 = 1, C, = 0 as usual, the term x = 0 coming from the pole at 
= 0). X converges everywhere. Hence fort < 7, expression (17) becomes 


— 1 

isina > (—1)"*7——S da sa(t(—1 "et OFEX; — (20) 
ae lt-- 3 
n=1 ” 


Gg 











446 M. J. LIGHTHILL 


and so when +(Q—R) encircles s = 7, in the positive direction (con- 
sidering (7, —@) as a right-handed system of axes, so that 6 < 0 comes 
first in an encirclement of (7,, 0) in the positive sense starting from + > 7,) 
it increases by —2X. Hence Q ‘has a period of +2X at +r = 7,’ (ie. 
increases by this quantity when the point is encircled in the positive 
sense), according as it has been continued from its starting-point via 
6s 0. 

We also have, for 7 <7,, 9S 0, with the positive integers the only 
poles to the right of the contour, 





ot 
ale ] T 1 — 
P = 27sin «a ——. : a a ob, (rer Fer 81-10) dv. (21) 
9 : v 
a7 $1N 7rv Vv 
—oi 


When 7 > 7, this is the sum of the residues at the poles to the left of the 
contour; these are all double. Hence, 


oe . d ] ae ; 
P = isina > (—1)"5 | (vtm)(1——) yp, (rein 
dv v 
n=0 — 
ba — d ] - 
= 7sINna > — | +m)(1—<)y,(rjevore a 
dy v ; 
n=0 v=—n 
PO ee l : i an 
+-7sin « > (Fim) 1+ -](—nC, (7) er, (22) 
n=0 . 
and the second term is +X. Hence, when P encircles 7 = 7, in the 
positive direction starting from 7 <7, (so that @ >0 comes first) it 


increases by 2X. 

Thus, if continued via 6 > 0 (right-hand sheet), P+ @Q has zero period 
at + = 7,, 6 = 0: in fact there is no singularity there at all since, by (21) 
and (18), P+@Q differs by a function regular at 7 = 7,, 6 = 0 from 


7a J 
sin « 








ed +4 | ne I ob, (r)eria—1-19) dy, (23) 


Qi cos7v = sin av 


— cot 
which converges uniformly in the neighbourhood of (7,,0) since the term 
in square brackets is 27ie-‘*” cosec2mv. Thus there is no singularity of 
P+ Q (except B) on the right-hand sheet: on the left-hand sheet, how- 
ever, there is a point S where P+ Q has period 4X. If a cut is made from 
this point to infinity P+@Q is one-valued on the cut Riemann surface, 
taking values as shown in Fig. 2, in the different regions. As fb = 3(P+4) 
encircles the point S (which corresponds of course to the point at infinity 
in the physical plane) it increases by 3(4X). By (19) this is not zero for 
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(con- | 4j] 7,8: hence the stream-function is many-valued in the physical plane 
omes | and the method (when f,,(7,) is taken as e-”*) is useless in this example. 
7 This is due to the presence of circulation. 
(i.e. 
sitive 
it via 
only 
(21 
Fic. 2 
of the ‘ 
4, The correct value of /,,(7,) 

To overcome this difficulty I thought of choosing a different f,,(7,) more 
symmetrical to %,(7). The reasons for this are far clearer from the 
method of ref. (2) which I evolved subsequently. _,,(7,) was first tried: 
this got rid of part of the period 3(4X) and pointed clearly the way to the 
orrect value, whose raison d’étre is demonstrated far more generally as 
well as convincingly in ref. 2 

] ) j 
: s_,(7,) +27, b_,,(7 
” falt,) = Cate Sana (24) 
(22 l—n 
where ys_,, denotes the derivative of %_,,. The work is here displayed with 
n the J this latter value used. 
st) it Consider, for 6 S$ 0, the integral 
a ] - (v—1)(—v—1)! (v—3)! ,, i 
eriod — ~ =~ F(v—}4, —4; $; sin’a)b,(7) x 
y (21 ” 4 (—2): 
~ \L9 , (2 
- wb, ( 7 ) T aT wb (7 ) ei? a—6) dv. (25) 
; l—v 
(23 for 7 <7, this is equal to the sum of the residues of the integrand at its 
poles on the right (these are double poles at the non-negative integers), i.e. to 
. term te n d | al? 1)(—v—1)! (vp—3)! 7 *¢ 
= i. = (~1) |» -n)? ; ) 2) F(v—4, —3; 4; sin®x)p,(r) x 
ity 0! n=0 - (—3): 
1 ¢ ' 
how- _ b(t) +27, f_ (71) _ 6 
J w i € —_ 
» from = 
f v=n 
iriace 
ak (—1)"—'(n—1)(n—3)! _— ~ 
PLO - ) (Fiz)(—1) “— F(n—4, —4; 4; sin®a)y,,(7)nC,, x 
. —_ n! L)! 2 Ne 
ifinity : , 
b (7) +27, u,,(7 
ro fot ye Unl7) 1Yin( ) pind (26) 
—n 
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which we write as Q+7', where 7’ converges everywhere and 7' + 0 as 
the Mach number at infinity tends to zero. Also, as this happens, the 


integral (25) tends to the low-speed integral (8), since y,(7) ~ 7?”, 
bola) +2nWoln) Pte (—br eA) yan 
l—v l—yv saps ve 


and (7/7,)!e-*”9 = [(q?/q2,)/(1/q2,) ve”? = (qe-)” = Cr. But (8) is the 
low-speed Q. Hence the Q here defined tends to the low-speed Q as the 
Mach number at infinity tends to zero and so is a reasonable extension 
thereof. 

When 7 >7,, the integral (25) is equal to the sum of its residues at 
4—n(n > 0) and —n (n D> 2), i.e. to 


iy (Foy : Mn 2): F(—n, —34; $; sin?a)ypy_,,(7) > 

| ai (— 2)! ble : 
: us, 4(7,) +27, ub, 4(7) 
“~~ 1 


n--s 


et(n—3)0 


Se ee 
> n ‘ts n- 2): p 


(—n—4, —4; 4; sin?«)(—nC,, u,,(7)) > 


w.\t9 P 
y b,,(73) aT} Yn(71) (__yyngind (28) 
n+1 

which we write as +W-+R, where R converges everywhere. Thus the 
continuation of Q for 7 >7, via 6S 0 is R+(T+W), and its continua- 
tion right round a cut from t = 7,, 0 = 0 to rT = 75, 6 = Ois 2R—Q; and 
continuation a second time round gives again Q. 

Now W can be rewritten as 
wo (—n—4)(n—3)![(—4)! nr! . 
iS , 2) = -“—_- sin a+ O(n cos?"x) | ys4_,,(7) > 

hot n!(—4)! (n—34)! 


Hy —4(7) +274 Yn4(71) 


n-+4 


wr 
“~N 


= ising > (4—n)—,_,,(7)(%, (7) +274 Wy_4(71) )e@-” (29) 
n=0 


plus a function analytic for 7 >7,. This series is 


cot 
= sina {| a7 u¢l(r Le aa ; 
hee ee 7 =5 (7) cy y(7) +27, yb_,(7,) er =i7- dy +X (30) 
2a7t_=—J ~CcOosmy op 
—< i 


according as 6 < 0, where 


X =7sina > C, u(r) (us,,(7,) +27, oi, (7,))e@”?, 
n=0 
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which converges everywhere. For + < 7,, therefore, (29) becomes 


x 


+sin »| > (n+3)7 vb, 44(7)(¢_, (74) +27, bi ,-4(74)) (File e— 
n=1 
b 3 mes, (7)C,, (eb, (71) +27, bi(rs))e-9] FX =— ZTY, say, (32) 
n=0 


where Y=asina > C,%,(r)(b,(7y) +27, hn(74)) (e+e), (33) 
n=0 
which is purely real. Thus W, when continued round + = 7,, 6 = 0 in 
the positive direction (6 < 0 first), decreases by 2Y: hence the period of 
Q= R+(7T+W) is +2Y according as it has been continued from its 
starting-point via 6 S 0. 
Consider next the integral 


: (15) p(x) Pooltadt ra Pola) erste dy, (34) 
V, 


27 J sinzv l—yv 


which when +t < 7, is equal to 
i — l b_.(7,) +27, b'_ (7 
~isina 5 (vn) }y,(7) 71) _— 1) 6-0 — 
Lut } l—y 
n 0 v nt 


- . = | l b (7) +27, vi, (7 , > 
~isinan S (J in)(1— }yp,,(r) nC, ¥nl7a) Pal 1) g~ind — P+(Y—X), 


n) l—n 


kod 
n 0 


(35) 
with X, Y as in (31), (33). The P here defined reduces to the low-speed P 
as the Mach number at infinity tends to zero, by arguments similar to 
those for @. When 7 > 7, the integral (34) is 


; . b_,(7,) +27, 0'_,(r sailed 
isin x ~ a ! n)(1- “) gh,(r) © y(71) nt 1p_,( 1) ,-iv8 


food «CLV j 
n 0 v ” 


bt ee l os (7) +27, oi,(7,) , — 
t+isina S (+ im)( + \( ne us,(7)) n 1) 71 Wy ) pind — HX say. 
n l+n : 
n=0 

(36) 
Hence P continues for r > 7, into H+Y, so that in encircling + = 7, 
’= 0 in the positive direction (6 > 0 first) it increases by 2Y. As in §3, 
therefore, P+-Q is one-valued on the Riemann surface cut to prevent 
encirclement of S, is regular except at S and B, and has period 4Y at S. 
¥, however, being the imaginary part of P+Q, has no period at S since 
Y is real, and so is one-valued on the uncut Riemann surface. The 
‘olution for ys is therefore quite possible physically; and an exact flow of 
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a compressible fluid past a contour with circulation has been mathe. 
matically specified. 

The series for % in different parts of the Riemann surface are indicated 
in Fig. 3. The value when 7, <7 <7, is obtainable by splitting W as 


$ S(H+R I(H+R 


+7+W) 


T=7, -T-W) 











Fra. 3. 


in (29) into a term continuable round 7 = 7, and one convergent for 
T >7>, and using ¢ = 3(P+R+FT+FW). 

Strictly, an investigation into the one-valuedness of x, y as well as that 
of ys is necessary, since x, y may have constant periods though ¢ has none. 
Further work is omitted here, however, since the required investigation 
is given in ref. (2). It is clear that the method given there is greatly 
preferable to the present one, as well for the reasons given in § 1 as for 
the fact that it gives a single series for % convergent throughout 7 < ,, 
which will be more useful than those of the present paper especially around 
T = 7, and 7, (series of the type of (47) of Part I will also be convenient 
in this region): this single series is also such that it does not involve 
relations in the low-speed hodograph plane, which our §2 shows to be 
rather involved in even a simple example. 
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THE METHOD OF CHARACTERISTICS FOR PROBLEMS 
OF COMPRESSIBLE FLOW INVOLVING TWO 
INDEPENDENT VARIABLES 


By R. E. MEYER 
(Department of Mathematics, The University, Manchester) 


PART II. INTEGRATION ALONG A MACH LINE. THE RADIAL 
FOCUSING EFFECT IN AXIALLY SYMMETRICAL FLOW 


[Rec ived 13 December 1947] 


L 


SUMMARY 

The paper deals with the growth and decay of disturbances along Mach lines in 
sentropic, irrotational, steady, two-dimensional or axially symmetrical, supersonic 
flow; in particular, the distribution of disturbances is investigated along a Mach 
ine in axially symmetrical flow on which the velocity is constant. As an example, 
the field of flow in the entry of a contractor of circular cross-section is calculated 
from the focusing laws, and the analytical expressions are compared with the results 
f the numerical methods of Massau and Tupper. 

The disturbance generated by the diffuser entry leads to a singularity of the flow 
pattern on the axis, the nature of which is investigated within the framework of 


inear theory. 


I, Integration along a Mach line 
TuE physical significance of the characteristics lies in their role as carriers of 
small disturbances. The whole pattern of a shock-free flow, for example, in 
aduct with uniform initial flow, can be interpreted as the cumulative effect 
of small disturbances generated by the walls, and the idea therefore 
suggests itself to study their growth or decay along the characteristics. 

The key to this investigation is found in the characteristic equation 
lla) of Part I,t as far as the Mach lines are concerned. If we employ 
general orthogonal coordinates based on one of the families of Mach 
lines, this equation is an identity, and hence the derivative with respect 
to 8 of the left-hand side also vanishes. The new equation thus arrived 
at governs the distribution along the Mach lines of the normal derivatives 
of the velocity components and the variables of state. 

If we restrict ourselves, for the sake of simplicity, to the case of an 
isentropic, irrotational, steady flow of a perfect gas, the new equation may 

t Cf. also ref. (1), § II, 11. 

t Ref. (2). Since many results of the §§ II, III, IV, and VI of Part I are used in the 


present report, the notation of Part I is followed throughout, and the equations are 
humbered through so that all numbers below (40) refer to the formulae of Part 1. 
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be written in the form 


dg 
pee SE eae +B = 0 4]) 
a ee a 
(the details of the calculation are given in the appendix (i)). 
The normal derivatives of the velocity components are related to the 


curvature of the normals by 


—] 
evg/hgeB = F(x) —*—v, kg, 
vg 1g OB (x) a 
4] (42) 
Ov, hg ep == vgkg— > (vg/4) F(a). 
dpe. 
(41) can also be written in the form 
d*h oh, \ dh 
Atte HE Pe nei, 3 nos Se he = ( : 
12 doe ( (2) aie (alg = 0, (88) 


and this equation governs the divergence or convergence of the Mach lines, 

The functions A, B, F depend only on the state of motion on, and the 
shape of, the Mach lines along which the equations hold; (41) and (42) 
therefore permit us to calculate the ‘propagation’ along a Mach line of 
a discontinuity of the normal derivatives of the velocity components as 
soon as the shape of, and the state of motion on, this Mach line are known. 
Unfortunately, A and B are of a very complicated form (cf. appendix (i)), 
and the equations (41) to (43) have so far been investigated only for Mach 
lines on which the state of motion is constant and which are therefore 
straight lines. In two-dimensional flow such a Mach line must belong to 
a region of uniform flow or to a simple wave, and an investigation of the 
focusing effect on it cannot contribute much to the known facts about 
these flow patterns. But for axially symmetrical flow we are led toa 
theory of the radial focusing effect which is the fundamental feature 
which distinguishes it from plane, shock-free, supersonic flow. 


II. The focusing of disturbances in axially symmetrical flow 
II.1. On a Mach line along which the state of motion is constant we have 
vg Ga = Ov,/0a = k, = 0, (44) 
and if we wish to include the axis in the region under investigation, we 


must also have R=vjr=0 (44a) 


_ 


(v, denotes the radial velocity component). Equation (41) then reduces to 


S (gseg) + ON, heey = 0 (45) 


Ox r 
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(cf. appendix (ii)), and since h,dasine = dr we have 
hgxg = mr-, (46) 


where m is a constant. By virtue of (3), a second integration gives 


C-- Tf" ° . 
hg = ; , with 1l+c = sine/2m, (47) 
where the constant of integration has been chosen so that hg = 1 when 
r=1. Finally, 
sine 
Kg = ——r7*(c+rt)-! (r #0), (48) 


and since F(«) = 0, owing to (44), equation (42) becomes 


ov y+l vB Ug aa 


= Ughgkg. (49) 


a” 4. e 

The rates of change of the velocity components per unit length in the 
direction normal to the Mach line are given by @v,/hgéB and dvg/hg ep, 
respectively. Both are proportional to 

r-*(c+rt)-1, 
The changes of the velocity components along the normal, from the Mach 
line to a neighbouring one, are given by (@v,/@8)dB8 and (évg/eB) dB respec- 
tively, and they are proportional to r-}. 

II.2. Equations (47) and (48) give hg and xg all along the Mach line 
as soon as either is known at any one point of it. If «g, the curvature of 
the normal, is discontinuous at some point of the Mach line, then xg is 
discontinuous at every point of it, and (48) determines the variation of 
the magnitude of this discontinuity. 

The curvature of the streamlines is discontinuous, too, where they cross 
a Mach line with such a disturbance. By virtue of Bernoulli’s equation 
we can write (12c) in the form 

1 ow 
“= ap hon’ 
replacing « and 8 by s and n when they refer to a coordinate system based 
on the streamlines. But stream and Mach coordinates are related by 
+thgd8 = h,dssinu+h, dncos p, 
where the upper sign corresponds, as always, to the upper sign in (11), 
and hence, by (44), 


cos ow 
a 
By (21), (11), (10), (44) and the isentropic relation 
po ay-d, (50) 
this becomes Kg + {2/(y+1)}cos*u sin p kg. (51) 
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Equations (51) and (47) allow us to draw some conclusions as regards 
the maximum disturbance which a wall may be permitted to introduce 
into a uniform, axially symmetrical flow without producing a shock-wave, 
Let us consider a straight pipe of radius r = 1, with uniform flow upstream 
of a certain point x = 0, where its cross-section begins to change in such 
a way that the radius of curvature of the meridian suddenly assumes a 
finite value R,,. The discontinuity of the streamline curvature persists 
along the ‘leading’ Mach line through the point x = 0, r = 1 which 
separates the uniform flow upstream from the non-uniform flow down- 
stream. It is the last Mach line on which the conditions (44) hold. On 
it, 6 = 0 and, by (16), « = -y and hence, by (49) and (48), 


R, = Un, = Neo+rirs, 


8 


, , +1)M4 
where N = (y+1)/sin*p cos?u = y 1) Mo >0 
Mi-1 
and M, denotes the Mach number of the uniform flow upstream. Thus, 
l+c= R,,/N, 
and by (47 ; 
and by (47), hg = 1—N(1—r')/R,,. (52) 


If R,, is negative, that is, if the disturbance is an expansion, hg is positive 
all along the leading Mach line inside the duct, but if it is a compression, 
hg > 0 only when r! > 1—R,,/N. When hg = 0 the leading Mach line 
meets an envelope of its family, and the presence of an envelope indicates 
the occurrence of a shock-wave. The inequality 


Ry > (vy +1)M4/(¢R—1) (53) 


is therefore a necessary condition for shock-free flow in the entry of a 
supersonic contractor of circular cross-section.t In the case of a given 
contractor, with small initial curvature, (53) determines two limiting 
Mach numbers between which shock-free flow is possible. The Mach 
number which admits the highest initial curvature of a contractor is 
M, = v2; and if the initial radius of curvature is smaller than 4(y+1) 
times the radius of the pipe, a shock-free flow is not possible for any Mach 
number. 

A similar condition can be deduced for two-dimensional flow. Instead 
of (46) we have hgxg = B = const., and if ¢ denotes the length on the 
leading Mach line measured from the wall, and hg = 1 at the wall, 


R, = +4N(t+B)sin p. 


t This condition has also been deduced by K. Friedrichs. 
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METHOD OF CHARACTERISTICS 
For a symmetrical channel of initial width 26, 
he ; 1-3N(b—y) Ry (54) 


on the leading Mach lines, if y denotes the distance from the axis. The 
condition for shock-free flow becomes 


R,, 4(y+1)M4b/(M2—1). (55) 


The Mach number which admits the highest initial curvature of a two- 
dimensional contractor is again, M, = v2; and if the initial radius of 
curvature is smaller than (y+ 1) times the initial width of the symmetrical 
channel, a shock-free flow is not possible for any Mach number. 
Equations (53) and (55) are, of course, necessary but not sufficient 
conditions for shock-free flow. 
A curious case occurs when R,, = N so that hg vanishes just on the 


axis. 


II.3. The question now arises whether similar results can be derived 
for disturbances of higher orders, that is, discontinuities of the curvatures, 
and normal velocity derivatives, of higher orders. Let us consider a Mach 


line on which 


CUp/Ca eee ae lip ep x= VY, 
CU ,/ Ca eee o ly /OB"da = 0, vas 
(56 
c=... = Oe /O0" = R= ... = ROP” = 0, 
x x i 
Kg mad om leg op*- —_ 0. 


Without loss of generality, we can choose the parameter, «, so that h, = 1 
on this Mach line. If we now differentiate (11 @) (m+ 1) times with respect 
to 8 and make use of (56), (11), (10), (4), and (50), we find that 


E™(hgxg)/EB™, o™+1yg/6B™+1, amtly /eB™+1 varyas rt (57) 
on the Mach line. If m 0, hg is given by (47), but if m > 1, 
eo. _ S-th Jome-t _. 
he l, Che = wa. Se hg ep == §, in 
(58) 


e™he/eB™ oc (r'—1) 
(the constant of integration is chosen so that @”"hg/é8" = 0 when r = 1). 
The equations (46) to (48) and (57), (58) are a generalization of a known 
result of linear theory (ref. 3). If we differentiate the equation for the 
linearized potential, 


tw 


] 
bnt+—¢,-—0%-p=0, of 


r 


= M?-1, 


(the suffixes denote partial derivatives) with respect to 2 and eliminate 
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¢,,. by means of the identities 
dd, = Perr di+byry dr, 
dd, fi Pre dx+- bir dr 
we alrive at 


(dx? — x? dr*)b,r7+db,, dr—d¢,, dx + : Pr dr? = 0. (59) 
r 


But on a Mach line where the state of motion is constant, we have 


dz+adr = 0, 
dd, = bre dx+y, dr = 0, 
- l d ' 
and (59) reduces to aj—+2—Ij¢,, = 0, 
r dr 
that is, Ging & #*, (60) 


provided M, #1. A similar result can obviously be established for the 
derivatives of any order. 

In a certain sense there is full agreement between the linear and the 
non-linear theory. If we interpret ¢,.. as the rate of change in axial 
velocity from one Mach line to another, the radial focusing law is the same 
in both theories; and this law does not depend on the order of the 
disturbance, provided all those of lower order vanish. The main improve- 
ment which the non-linear theory brings about is the distinction between 
the rate of change from one Mach line to another and the rate of change 
per unit length. It is an unexpected result that the error to which the 
linear treatment may lead if these two notions are confused increases in 
importance, not as we approach the axis, but as we go away from it. 

If we consider the entry of a circular duct with shock-free flow our 
results may be summed up as follows. If we calculate the velocity com- 
ponents on a Mach line which belongs to the same family as the leading 
Mach line we obtain the same result from both theories, provided we take 
the ‘second’ Mach line near enough to the leading one for the velocity 
components to be given by the linear terms in their Taylor series with 
respect to 8.t This is true for any ‘smoothness’ of the wall. But the 
simple Mach line pattern of the linear theory is distorted in the general 
theory, and the distortion depends on the smoothness of the wall. 


III. The entry of a supersonic contractor of circular cross-section 


III.1. The theory of radial focusing enables us to explain the surprising 
results obtained by Tupper in his paper on a supersonic contractor (ref. 4) 


+ It will be seen from the results of § IV that this assertion cannot be extended to 4 
certain portion of the second Mach line near the axis. 
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to which this investigation owes a good deal of stimulation. Tupper’s 
work was based on the linearized equations of motion, and as far as the 
linear theory goes the explanation of the flow pattern in the entry of an 
axially symmetrical duct will be completed by the results of § IV. 

In order to investigate how far the focusing laws given in § II.1 can 
explain the flow pattern which the non-linear theory predicts for the 
entry of such a duct, a characteristic very near to the entry of Tupper’s 
contractor was recomputed, using the full, non-linear equations of motion 
as described in § VI and appendix (ii) of Part I (ref. 2). The result is 


shown in Table 1. The initial Mach number (M, = 1-9) and the shape 


TABLE 1 


Second Characteristic of Contractor 














0 Vx Vp (70) (71) 
5 2 0437 | 1°5517 —*o006gI —*006878 —*OO7II5 
0460 | 1°5815 —*00727(5)| —*007244 "007512 
0455 I°5512 —*00772 —*007675(5) "007983 
7 2 521  1°5509 —'00524 —*008195 "008556 
7 562 | 1°5504(5) | —*00585 —*0055 39 *009273 
2 615 | 1°5799 —*0097 2 —*009665 "010206 
7 650 | 1°5792 "01053 —*O10730 *O11491(5) 
{ 791 | 1°5780 "01245 —*O12407 "013428 
67 | 1°5762 —"O1524 —*O15105 *O16850 
55 —*OIII7Z | 1°5745 —°o170 —°01735 *019954 
+7 7 37 Is7is -"O215 -*02094 —*025795 
7 $595 | °7 165 1*5680 "0259 —°02417 —*033343 
7 $301 | °7 21¢ 1°5604 "0337 —*02207 —*054686 
404 7 260 | I°5501 "0403 +*10043 —*IOgIoo 
524 °7 —"0258 1°5332 —*0390 
TO0§ | °702 1°4797 ° 
Wa 0°0312740x? —0:00064855x°. 
I am } SS ’ 2°5156794, as 0°6968641, Uo wg = 1°5860893. 


of the wall, (62), are the same as in Tupper’s paper, y = 1-4 and the 
second’ Mach line was chosen to meet the wall at 2 = 0-07 so that it 
starts very nearly, at least, at the same place as the characteristic m = 2 
in Tupper’s case. The table gives the coordinates, x, and rz, of the points 
which were computed together with the values, at those points, of the 
squares of the velocity magnitude, w, and the speed of sound, a, of the 
ingle, 6, between the velocity direction and the axis, and of the radial and 
axial velocity components, v, and v,. All velocities are understood as 
nultiples of the critical speed of sound, a,. 

According to the linear theory, the second Mach line is straight and 
parallel to the leading one. The present computation shows that it is, in 
lact, nearly straight; but the distance between the two Mach lines decreases 
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appreciably as the radius decreases. The leading Mach line meets the axis 
at x = ,/(M%—1) = 1-61555, and the table therefore shows that this 
distance (measured in the axial direction) falls from 0-0700 at the wall 
to 0-0227 at the axis. Since it is the radial velocity component which 
shows the most unexpected behaviour, according to Tupper’s graphs 
(ref. 3), it has once more been plotted in Fig. 1 according to the results 


Ved 
(71) 


--04 





all \ 























Fic. 1. Distribution of the radial velocity component on the second Mach line 
according to (i) the computation by Massau’s method, (ii) equation (70), and 
(iii) equation (71). 


of the present computation. A comparison shows that the rapid variation 
and the peak in the distribution of v, near the axis are even more notice- 
able than in Tupper’s results. 

We now proceed to compare the results of our computation with those 
obtained for the radial focusing effect. If we choose our system of general 
orthogonal coordinates so that 8 = 0 on the leading Mach line, the 
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equation of the second Mach line will be 
Bp = const. = B, < 1, 
and the radial velocity component on it is given by the Taylor series 
; . Q ; R 2/(2,, /AR2 ! ‘4 

Ure B,(cv, CP)g=o7 $P3(6°v,/68 )p=ot ++ (61) 
since v, vanishes on the leading Mach line. The variation of év,/é8 on the 
leading Mach line is determined by the focusing law for first-order 
disturbances, and the equations (49) and (46) show that 

(ov, ©B)g aq r—t, 

Clearly, a part only of the distribution of v, shown in Fig. 1 can be approxi- 
mated by the first term of the series (61). However, we must not forget 
that equations (46) to (49) only describe the focusing of the first-order 
disturbances, whereas the wall, the shape of which is given by the equation 

l—r = da?— 0-0045044824-+ 0-00021643026 (62) 


(cf. ref. 4), introduces also discontinuities of the third- and fifth-order 
curvatures of the streamlines at the point 2 = 0. 

The calculation of the focusing laws for these higher order curvatures, 
in a case where the lower order ones do not vanish, would be a tedious 
task, and we shall try to avoid it. The description of the first step of the 
Massau process given in § VI.1 of Part I (ref. 2) shows that only the 
following data about the wall enter into the computation of the second 
characteristic: the values of x, r, and @ at the points where the leading 
and second Mach line meet the wall. With our choice of the (z,7r)- 
coordinate system, and with an equation of the form 


l—r = Ax?+ Br3+ Cxt+ Drd+... (63) 


for the shape of the wall (so that r = 1, 0 = 0 at x = 0) we can start the 
second characteristic with any proposed values of x, r, @ by an appro- 
priate choice of two of the coefficients, while all the others are put equal 
to zero. For our purpose, a wall of the shape (62) is equivalent to one of 
the shape rae a i ene ; 

l—r 0-03127402?— 0-000648552°. (64) 
In fact, even if a contractor was designed in such a way as to avoid any 
discontinuity of the first-order curvature of the wall at the entry, we 
should still be justified in computing its second characteristic as if the shape 
of the wall were of the form (63) with C = D =... = 0. 


III.2. We therefore propose to evaluate the first two terms of the series 
61) and to compare the result with that of the computation. The radial 


velocity is 
, v UV, SIN €+- Ug COS € 
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ov ov, . ov F Ce 
and hence Se ‘tn. a hee cos e+ (v, COS e—Ug sin €) — 
6p op op ep 


= —— whgKp cos? 
y+l oF 
on the leading Mach line, by virtue of (49), (4), (16), and since v, = weosy, 
Taking account of (46) and (47) we therefore have 
weos*usine _, - 
——_— 7, #, (65 
(y-+1)(1+e) | 


where the suffix 1 indicates that we mean the radius of the point on the 


(ov, €B)g—o _— 


leading Mach line which lies on the same normal as the point on the second 
Mach line at which we propose to calculate v,. 

In order to find the focusing law for the second-order disturbance we 
have to differentiate (11a) twice with respect to P taking account of 
(44). We do not now assume, as we did in § II.3, that the first-order 
disturbance vanishes (i.e. xg, etc. are now + 0), and a somewhat lengthy 
transformation of the same type as that given in the appendix (i) is needed 
in order to deduce an equation for the variation of @(hgxg)/éB with « along 
the leading Mach line.+ It is found to have the form 

d (hg kg) sine O(hg kg) 


sine 3+ _ dy—l1 ela ~ 
“ Esies Lo +. — 1 $ tan ¢ -+-——_- 60 h 24 
h, da ep 2r . OB ' 4r "a y+l1 €}( Bp) 


4 = jtan e+. (4 Pe... 11)eo clip = 0. (66) 
Sr* | y+1 ; 
In contrast to (45), this equation is inhomogeneous, which means that 
a first-order disturbance at the entry of the duct always generates a 
second-order disturbance even if the equation of the wall, as, for example, 
(62), shows no discontinuity of the second-order curvature of the wall. The 
algebra which leads to (66) makes it appear highly probable that the 
equations for the variation of the disturbances of any order will be found 
to be inhomogeneous whenever a disturbance of lower order is present. 
We must expect a disturbance of any order to generate disturbances of all 
higher orders. 
If we introduce (46) and (47) into (66) and note that dr = h,dasine 
on the leading Mach line, we can write (66) in the form 


d é ; i . 
al gators) = m?(A,r-2—B, r-), 


74 — 15 ay —- 4 
“Y** cot e—tane}, B, = tane —*—"— cote. 
y+l 
+ This algebra has been omitted through lack of space. A detailed account of it will 
be found in the appendix of A.R.C. 10649. 
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This integrates to 
af 
op 


In order to find @hg/éP itself we note that, by (3) and (44), 
Cxg he hg @ [ eh hg @ 
h 2 *B - . (h . Kp) = B > B = B er log hg, 
“Op =h, op ‘ h,, OB \hg ea h, @a0B ~ 


a 
and therefore, again by virtue of (3), 


he Che h he 0 oh 
75 ¢ : op 


op OB hgh, ea h, ex h, @x 


hgxg) = m*(2B, r-1— A, r-#+ 6r-4). (67) 


It follows that 


] l f 
; 8 ~ —(hgxg) dr 
op sine J op * * 
1 
2m? 1 iN } ve] 
—[A,(r+—1)+ B, log, r—8(1—r')], (68) 
sine 


if the lower limit is chosen so that éhg/eB = 0 when r = 
From (67) it can be deduced (cf. appendix (iii)) that 


ev 2m* 
r - 1cos2 - $1 Set p 
(- ed = ——-o cos*e[ B, r-1— A, r-?+-dr-], (69) 
B*] 2-9 ) 
; ( 3y+7 
with A, ; (3 tane 11 cote), B, = —tane— we - cote, 
i 


and combining this with (65) and (61) we have 


— 
Ww COS“e SIN € , —— _ 
Vig = B.{{4(1+c¢)+f,8 sin e}ry?+ 


" 4(y+1)(1+e)? 

| +-Byrz1sin e(B,—A,rz*)]+0(68), (70) 
provided 1, 0 (since 6 0, ~= +e on the leading Mach line; 
m= 4sine/(1+c)). The constants of integration, c and 5, are found by 
comparing (65) and (69) with (64), at r, = 1.¢ B, can then be evaluated, 
by the help of (47) and (68), from the distance between the leading and 
the second Mach line, measured along the normal through the point of 
intersection of the second Mach line and the wall. The values of the 
constants which correspond to the initial data of the computation are thus 

found to be 

c = 03341401, 4 0-9108727, PB, = 0-0371515. 

In order to compare the distribution of v, on the second Mach line which 
(70) predicts, with that obtained from the computation, the values of 7, 
which correspond to the point 29, 7, given in the table were computed with 


+ Cf. last footnote. 
Hh 
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due regard to the slope and curvature of the normals; v,. could then be 
evaluated from (70) as a function of r,, and the result is found in the last 
column but one of Table 1. For further comparison, v,. was evaluated 
from the equation a 
weos’esine , _ 
‘4 = s—= — Bar; *, (71) 
(y+1)(1+-e) 

which is obtained from (61) and (65) if account is taken only of the first- 
order disturbance at the entry of the contractor, and this result has been 
entered into the last column of the table. 

The distributions of v, as found from (70) and (71) are also plotted in 
Fig. 1, and it is seen that the rough approximation obtained from (71) is 
much improved by taking account of the focusing effect for the second- 
order disturbance. The most striking feature of the curve (70) is that it 
shows indeed a peak, owing to the higher negative powers of 7, in the 
‘small’ terms of the formula. 


III.3. Our comparison shows that over a considerable part of the entry 
of an axially symmetrical duct, the field of flow can be approximated very 
satisfactorily by analytical expressions derived from the focusing laws 
(46), (47) and (67), (68) for the first- and second-order disturbances. 
However, as the radius decreases, a systematic discrepancy appears, and 
the agreement becomes unsatisfactory for small values of r. 

Before we accept this result it may be useful to check the accuracy 

of the numerical method. In any computation small random errors are 
introduced at every step, and the question arises whether they tend to 
cancel out in our case. In order to investigate this question we shall 
imagine that a small error is introduced into the value of one of the 
variables at some point of the second characteristic, and that the com- 
putation is then continued once on the basis of the correct values at that 
point, and a second time on the basis of the erroneous values, assuming 
in both cases that no further error occurs. Instead of comparing these 
two computations we can compare the second characteristics of two con- 
tractors which have initial radii roughly equal to the radius at which the 
error was introduced, and suitably chosen small differences in the initial 
Mach number, the initial radius, and the shape of the wall. We now apply 
the equation (71). There will be small differences in the values of the 
constants on its right-hand side, and it follows that the absolute difference 
in the values of v,. increases as rj + while the relative difference is inde- 
pendent of the radius. 

+ The argument can be followed up in more detail, including, for example, the modifica- 


tion of the position of the second Mach line which is caused by the introduction .of the 
error. Cf. Appendix B of A.R.C. 10649. 





Th 
focus' 
than 


comp 
crepa 
veloc 
focus 
the ¢ 
stren 
rate 
in th 
tive 
betw 
TI 
the 1 
the 
deri 
of tl 
of tl 
cale 
axis 
ana 
tion 
the 
In f 
poir 
vari 
rant 





en be 
e last 
uated 


(71) 


first- 
been 


ed in 
71) is 
cond- 
hat it 
n the 


entry 
| very 

laws 
inces. 


, and 


racy 
"Ss are 
nd to 
shall 
f the 
com- 
that 
ming 
these 
. con- 
h the 
nitial 
upply 
f the 
rence 
inde- 


\difica- 
of the 











METHOD OF CHARACTERISTICS 463 


The argument shows that small computational errors are, in fact, 
focused like disturbances, and a total error which is considerably bigger 
than the sum of the small random errors introduced at every step of the 
computation must therefore accumulate from them. A systematic dis- 
crepancy must exist between the computed distribution of the radial 
velocity on the second Mach line and the true distribution. However, the 
focusing of errors can cause a breakdown of the Massau method only in 
the cases where infinitesimal disturbances are actually focused to finite 
strength.+ In our case we cannot expect the errors to increase at a higher 
rate than the radial velocity itself, and since the random errors incurred 
in the present computation are of the order of 0-1 per cent. of the respec- 
tive values of v,, at most, their focusing cannot account for the discrepancy 
between the curves of Fig. 1. 

The improvement of the approximation which was brought about by 
the use of equation (70) instead of (71) suggests that the range over which 
the field of flow is represented satisfactorily by analytical expressions 
derived from the focusing laws can be extended farther by taking account 
of the third and higher order disturbances which are generated by those 
of the first and second order. However, our analytical method for the 
calculation of the second characteristic cannot be extended right to the 
axis. The normal which meets the leading Mach line on the axis constitutes 
anatural limit for the use of Taylor series of the type (61) for the calcula- 
tion of the dependent variables. Comparison of (71) and (70) shows how 
the convergence of the series deteriorates as this limit is approached. 
In fact, all our formulae indicate a singularity in the field of flow at the 
pont where the leading Mach line meets the axis, and the dependent 
variables are probably not analytic functions of 8 at all over any finite 
range, however small, on the normal through this point. 

It appears desirable, therefore, to investigate the singularity directly. 
A first step in this direction is described in the following section. 


IV. The singularity on the axis in linear theory 


IV.1. If the condition (53) is satisfied, we shall not expect that a dis- 
continuity of the velocity components themselves is caused by the focusing 
eflect on the leading Mach line, and the computation provides no con- 
trary evidence. Moreover, it has been shown in § II.3 that the agreement 
between the predictions of the linear and the non-linear theory improves 
as the axis is approached along the leading Mach line, and we shall there- 
lore try to find a first approximation for the field of flow near the singu- 
larity by the help of the linear theory. 


+ Similar conclusions appear to have been arrived at also by Tollmien (ref. 5). 
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The linearized potential which describes the flow pattern predicted by 
linear theory must be a solution of the equation 


] 2 
Pret be 0b x = 0, = My—1, (72) 


which vanishes upstream of, and on, the leading Mach line, and the second 
derivatives of which vary like r-? on the downstream side of this line and 
the first radial derivative of which vanishes all along the axis. The 
expression 

$ 2 


37 


| (x’+ar cos u)'5 du, (73) 
0 


with x’ = x—a, 6 = 1 when 2’+arcosu > 0, and § = 0 otherwise, satisfies 





these requirements. That it is a solution of (72) is known (ref. 3). The 
conditions for ¢ and ¢, are obviously satisfied, and the behaviour of 
(73) on the downstream side of the leading Mach line is found by 
introducing gq = (x’+-ar)/2ar and substituting sin?(u/2) = q sin?@, whence 


7 
Jo 1s 
$, =- | (a’ + ar cos u)*d du 


cos?6(1—qsin?@)-8d8+0 as q>0, 





7/2 


= —(2ar)-3 [ (l—qsin?@)-!6 d@ > (8ar)-? as gO. 
0 
This transformation also shows that the velocity components vary like 
r+, and their derivatives like r-!, on the lines where x’ /r is constant, that is, 
on the straight lines through the singular point 2’ = 0, r = 0. On the 
axis, downstream of this point, 
$, = Vz'.T (74) 
In order to find the behaviour of the solution on the ‘reflected’ Mach 


line which slopes downstream from the singular point, we substitute 
oS 
p = (ar—2’)/2ar and o = (7—u)/2, whence 


d > — (2ar)!, ¢,—> (2ar)}, dre ro as p>t0. (79) 


+ An indication of this result is to be found in ref. (3). 
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The behaviour of (73) near the axis, the reflected and the leading Mach 
lines, is best found by writing the potential as the product of r! and a 
function of 1/p, p, and q, respectively, and substituting in (72). This leads 
to a hypergeometric differential equation for each of the functions; in 
particular, if we write 


d= 2 (2aur)?K(p), 


> 
v 


4 


7 


(72) gives 


p(l—p)K"+(2p—1)K’—{K = 0. 
Near p = 0, the general solution is given by a series of the form 
K = a,+a,p+a, p*log p+ O(p*);T 
q, and a, are found from (75). Near the reflected Mach line the velocity 
derivatives therefore vary like 
i + 
r—tlog p.t 
The disturbance on the leading Mach line is reflected from the axis not 
as a disturbance but as a logarithmic singularity.§ 
From (75) and the corresponding value of ¢,, the series for ¢ valid near 
the reflected Mach line is found to be 


d (2ar)?(1—?p— p*log p+ O(p?)), (76) 


and similarly, the series near the leading Mach line and near the axis are 
found to be = F 
ad 1(2ar)? F'(4, 4; 3;q), (77) 
; ; : 4d 
d = 3(x’—ar)t F(t, —3; 1; 1/p), 


where F denotes the hypergeometric series. 


IV.2. In order to check whether the linearized potential, (73), really 
provides an approximation for the general, non-linear field of flow near 
the singularity, we must insert the expressions (76), (77), in the non-linear 
lifferential equation, 


9 
a* 
m2 2 9) . 2 2 | - 
(a2—@? )”., 20 OD ®.+(a @?)0 + 0. 0, 
2 2, Y—l),,2 22 
a a+; (wo —O;—O®>2), 
Ui. ref. (6), § 8 
This result has been confirmed, and extended, by Ward (ref. 7). 
§ With reference to Fig. 1, it may be noted that (76) does not imply a cusp of the velocity 


ition along the second Mach line; its point of intersection with the reflected Mach 


axis and the point where v, is a maximum. 
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466 
for the velocity potential ® = w,(x-+¢). It is seen that the non-linear 
terms are, in fact, small compared with the linear terms when 7 is small, 
provided p is not small. But, however small we may choose r, there is 
always a range of values of p for which this is not true, owing to the 
logarithmic singularity of the second derivatives of ¢. 

Moreover, if one tries to construct a closer approximation for the non- 
linear potential in the form of an expansion in powers of r with coefficients 
depending on p, with the first term given by (73), one finds that the 
expansion fails to converge if p is sufficiently small. The methods em- 
ployed in this investigation do not permit us to draw any conclusion as 
regards the flow pattern in the immediate neighbourhood of the reflected 
Mach line and downstream of it. 


V. Conclusions 


It is shown that the growth and decay of disturbances along the Mach 
lines is governed by an ordinary differential equation of the first order 
(eqn. (41)). The convergence and divergence of the Mach lines is governed 
by an ordinary, homogeneous, second-order differential equation (eqn. 
(43)). 

These equations are integrated for the case of a straight Mach line in 
axially symmetrical flow which occurs, for example, at the entry of a circular 
duct with uniform flow upstream. The discontinuities of the normal velocity 


1 


derivatives of any order vary like r~} along such a Mach line, provided 
that no discontinuities of the derivatives of lower order are present. This 
result is compared with the focusing law of linear theory. On the other 
hand, a disturbance of any order must be expected to generate disturbances 
of all higher orders. 

A disturbance of the first order is equivalent to a discontinuity of the 
curvature of the streamlines where they cross the Mach line in question. 


It follows from the focusing law that the initial curvature of the wall of 


a supersonic diffuser must be limited if the flow is to be shock-free (eqn. 
(53)). For a contractor with given initial curvature of the wall there is 
only a certain range of Mach numbers for which a shock-wave can be 
avoided. (A similar condition applies to the initial radius of curvature 
of a two-dimensional contractor (eqn. 55).) 

As an example, in § III, the entry of the contractor is calculated 
which was investigated by Tupper (ref. 4) with the help of a numerical 
method based on linear theory. The field of flow is computed by the 
Massau method, as described in Part I (ref. 2), and it is shown that a fair 
approximation can be obtained for the greatest part of the contractor 
entry by using analytical expressions which take account only of the 
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disturbances of the first and second order. Moreover, it is seen that the 
numerical method suffers from the disadvantage that small computational 
errors are also focused like disturbances. On the other hand, the field of 
application of the analytical method is limited, near the axis, by the 
normal which meets the leading Mach line on the axis; the convergence 
deteriorates as this limit is approached. 

The focusing law indicates a singularity of the field of flow at the point 
where the axis is met by the leading Mach line, which carries the 
disturbance. It is shown that, on linear theory, the field of flow is repre- 
sented by the linearized potential (73), which is evaluated in terms of 
hypergeometric functions. The velocity components are seen to vary like 
, and their derivatives like r~!, on all straight lines through the singular 
point. The discontinuity of the velocity derivatives on the leading Mach 
line is reflected from the axis as a logarithmic singularity. If the velocity 
components are plotted along any line which crosses the leading Mach 
line and its reflection these curves suffer a sudden change of slope at the 
former place and they have a vertical tangent (but no cusp) at the latter. 

Upstream of the reflected Mach line this solution provides an approxi- 
mation for the flow pattern of the general, non-linear theory near the 
singular point. However, the method fails to provide any information 
about the field of flow in the neighbourhood of the reflected Mach line. 

I should like to thank Professor Goldstein very sincerely for his en- 
couragement and helpful comment. 
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APPENDIX 

. d 2 dy, 

i) From (50) and (11), P = lve 
p Y l v2 


B 


nd eliminating p from the equation of continuity, (10), we have 


bo 


td c Cte 
| —ty Ketvek,+R 0 (78 
hea a 12 ] ’ h,eB xB Bal ’ ) 
wit} ay - vw 
th A (y+1)/(y—1). 
For a perfect gas, Bernoulli’s equation can be written 


vi+Avz constant 
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(by virtue of (11)); and hence 


v,dv, —Avgdvg, (79 
and (78) becomes 6v,/€B = hg F(a)—(v,/A)hgkg, (80) 
. VY; 2 v,\ ov ] 
} m = ( sae a cS 4. R): P 
oat Ma) Ve 7 l a x ek . (81) 


whence (42) follows. Moreover, Bernoulli’s equation and (11) and (79) may be used 
to eliminate p, p, a, and v, from ‘Spain which thus assumes the form 

If we now differentiate (82) with respect to B, (80) with respect to a, and eliminate 
*v,g/Ca0P and also eliminate @h,/éB by (3) and éx,/éB by (5), we arrive at the follow. 
ing equation: 





eK ‘ (viK F\ [2v, 2w? CU, 
+4 ad) +(§22-2) it a+ (or gaa 2( oni) 
Avg vg/ Lv, (y—1)v% P. h 
Qk, Vp Ov veg OR wu? ¢ I CUg 
4 Ste rte Pte) ve Oat 2 (Fe, en) 0, (83 
wv. v, h, ea vi hgop* : v, \h, €a h, ex 
with Q=2 ed~2 ly, 
2 ] : 
Now R = —(v, sine+ vg cose), 
= 
and ér/hgéB = cose, 
and hence, by virtue of (79), (80), and (4), 
cR 2 v.cose F : R 
— = <—ion-snel.” la a -(v, cos e—Avg sin €) — ~ COS €, (84) 
op y+ r rv, r 


which shows that (83) is, in fact, of the form (41). Equation (43) follows from 
(41) if we note that, by (3), 


a 





= haha BK} 3) thax 


(ii) Equation (45) is best obtained by taking account of (44) when differentiating 
(82) with respect to 8. Moreover, it follows from (44) that F(a) = 0, and hence 
that (49) holds on the leading Mach line and that 


v2 /CaeB —(v,/A)e(hgkg)/éx. (89) 
In order to eliminate éx,/é8 we note that, by (5), (44), and (3), 
éx,,/6B C(hgkg)/h, Ca; (86) 


finally, @R/eB can be eliminated by the help of (84) and (44) and we are left with 
COB € 


2v, O(hgkg)/ex- hg > eke 0, 


from which (45) follows, since R = 0 implies 


v,/v, = —tane (94) 
on the leading Mach line. 
(iii) It also follows from (87) that 
205 2p, 22), ‘ 2 
ont oc“ O-Vea cv c“€ 
_ St B x - B.: ° gg) 
+ = —“sine+—* cose- a ene sine) 4+(y cose—vgsine)—, | 
cB ope Bs Tape cose + 2c alae ep We sine) op 


on the leading Mach line. By differentiating (81) with respect to B and taking 
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account of (85), (49), (86), (84), (45), and repeatedly of (44) and (87), one finds that 





oF v, COSE j 
= la Kas 
cB (y t 1)Ar ane 
and hence, from (80) and since F 0, that 
C2u, v Olheke) Ve V,COSE 16 
e=_ — : PP heke  ... = ha 
op A op A \ “i e) (y+ 1)Ar oe 


on the leading Mach line. From (79) and (87), 

o [ ov 
é2v, /eB? = A- (tane—*), 
»/ ©} op\ ep) 


and if account is also taken of (4) and (49), (88) assumes the form 





Cn Qu Oheke) 4v, . COS € : v hk 
Sf on t..o08 € ——— —sin €(hg kg)? -+-——| vg sin e—-* cose EB | 
cB? y+ ] op Y ] , y+1 A r 


whence (69) follows by virtue of (87), (67), (46), and (47). 











ASSESSMENT OF ERRORS IN APPROXIMATE 
SOLUTIONS OF DIFFERENTIAL EQUATIONS 
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SUMMARY 


The term assessment is applied to any process which enables us to set rigid 
bounds to the error or to estimate its value. It is shown that upper and lower 
bounds can be assigned whenever the Green’s function of the problem is one-signed; 
this is true in many important problems. Another method is applicable to step-by- 
step solutions of ordinary differential equations, linear or non-linear, and depends 


on using the ‘index’ of the process of integration. Lastly, the error in a linear 
g : 


problem can be estimated when an approximation to the Green’s function is known. 


1. Introduction: Meaning of assessment 

THE term error is used to signify the difference between the value of an 
unknown as given by some process of approximation and its true value, 
When the differential equation is written as an expression equated to 
zero, the value of this expression for the approximate solution is called 
the residual. The word assessment is used to cover any process which 
enables us to estimate or delimit. There are two principal kinds of 
assessments of error: 

(a) the fixing of.rigid upper and lower bounds to the error; 

(b) the estimation, more or less closely, of the error. 

When it is possible to fix rigid bounds to the error these are usually rather 
widely separated and the assessment is correspondingly crude, but when 
the bounds are both rigid and close we have a most useful assessment. 
However, when we have not discovered how to set rigid bounds, or when 
these are insufficiently close, we must have recourse to some method of 
estimation. Such an estimation may yield results of ample accuracy for 
the purposes of applied mathematics even when lacking ideal precision. 
It is to be remarked that we exclude strict evaluation of the error from 
consideration since, if this could be done, we should have a method of 
exact solution. An essential requirement is that any method of assessment 
shall be applicable throughout the whole region of integration. 

We shall here confine attention to differential equations with boundary 
conditions which render the solution unique. No attempt at a general 
treatment will be made, but the following items will be discussed: 

1. Bounds to the errors of ordinary and partial differential problems 
having one-signed Green’s functions. 
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2. Estimation of the errors of the step-by-step solutions of ordinary 
differential equations, or sets of these, with one-point boundary 
conditions. 

3. Estimation of the errors of linear problems when an approximation 
to the Green’s function is known. 

Item 1 is concerned only with linear problems and may appear of very 
limited applicability; in fact it covers some very important problems such 
as Poisson’s equation with fixed boundary values of the unknown. The 
basis of the assessment is the value of the residual, and it is assumed that 
this can be found from the approximate solution. For approximations 
such as those provided by the methods of Rayleigh—Ritz, of Galerkin, or 
of the thickness parameter the residual will be given as an explicit func- 
tion, which is particularly convenient. Item 2 covers non-linear problems. 

The intention of the writer is to give a concise account of some useful 

methods and no attempt is made to enter into the minutiae of numerical 


work. 


2. Bounds to the errors in problems having one-signed Green’s 
functions 
2.1. Outline of the method 
We consider an ordinary or partial linear differential equation with 
boundary conditions of the linear and homogeneous type which render 
the solution unique. We suppose further that the Green’s function of the 
problem is one-signed. 
Let the differential equation be 
Ad = f, (2.1,1) 
where A is a linear ordinary or partial differential operator, ¢ is the 
unknown, and f is a given function of the independent variable or variables. 
Also let 
Ac = 1 (2.1,2) 
with the same boundary conditions as for d. We shall call o the basic 
solution, and it is one-signed on account of our assumption about the 
Green’s function. Suppose that ¢, is some approximation to ¢ which 
exactly satisfies the boundary conditions, and let 
Ad, —f = «, (2.1,3) 
so that ¢ is the residual corresponding to ¢, and is, in general, a function 
of the independent variables. Let 


€, = absolute maximum value of e in the region of integration, 


and 


absolute minimum value of « in the region of integration. 
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Then (e—e,) is everywhere negative or zero and (e—e,) everywhere positive 
or zero in the region. 

For definiteness let the Green’s function be everywhere positive in the 


region. Then by equations (2.1,1)—(2.1,3) 


A(¢—4¢,+4& oc) €,—€ = DY, 


and consequently —¢,+«,0 > 0, 

or d > ¢,—«19- 

Also A(@—¢,+€.0) = e—e <0 

and d < b4—€g 6. 

Finally, $,—§ 5 <b < b,—€gG. (2.1,4) 


If the Green’s function were everywhere negative we should have 


dp 


b,—€29 <6 < o,—E\G. (2.1,5) 
These inequalities call for two remarks: 

(a) Provided that ¢, and e, are small, a rough approximation to o will 
usually give sufficient information about the errors in ¢,. 

(b) A large but highly localized error in ¢,, or its derivatives, will yield 
numerically large values for one or both of €,, €, and will greatly 
widen the bounds in the inequalities. It is therefore important to 
avoid such large local errors when the present method of assessment 
is used. 

2.2. Errors in the basic solution 

Let o, be an approximation to o and 

99 
Ao,—1 = ». (2.2,1) 
Also let 7,, n2 be the absolute maximum and absolute minimum values, 
respectively, of 7. Then 
A(o—o,+7,¢) = 1-7 2 9 


and A(o—o,+72°) No—y <9. 


Hence, if the Green’s function is positive, 





Co , Co 9 ¢ 
Se <o< ——t_. (2.2,2) 
1+, 1+» 
but if the Green’s function is negative, 
oq - i oq (2 ) 3) 
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2.3. Cases where the boundary conditions are not homogeneous 

When the boundary conditions are linear but not homogeneous we can 
reduce the problem to one with homogeneous conditions as follows. Let 
8 be a convenient function which satisfies the boundary conditions. Then 


b= o—B (2.3,1) 
satisfies linear and homogeneous boundary conditions and 
Ay = f—AR, (2.3,2) 


where the function on the right-hand side of the equation is known. 


2.4. An important case where the Green’s function is one-signed 

By way of example we shall show that the Green’s function for Poisson’s 
equation is negative, the region of integration being the space interior to 
a closed surface B upon which the solution is to vanish. 


Let the differential equation to be solved be 
V4 = o, (2.4,1) 


where ys is everywhere positive within B and ¢ is zero on B. Then either 
¢is everywhere negative within B or it is positive at some point P within 
B. If it is positive at P it must be also positive everywhere within a closed 


surface C' surrounding P 


and vanish on C; moreover, C is entirely within 
B or coincides with it wholly or partly and % is therefore everywhere 
positive within C. Apply Green’s theorem to the space within C. 


[Lf (2) (2) + (2) +904} aeayas = — [ f sas =o, 


cy C2 OV 


(2.4.2) 
where €¢/év is the rate of change of ¢ along the inward drawn normal 
to C and dS is the element of the surface of C. But according to our 
hypothesis the integrand is everywhere positive within C and the integral 
cannot vanish. Accordingly ¢ cannot be positive anywhere within B, 
and, since % is an arbitrary positive function, it follows that the Green’s 
function is always negative. The same theorem is true in two dimensions. 

This proposition was applied by the writer some years ago to delimit 
the errors in approximate solutions of special problems in the theory of 
elasticity (1,2). In place of the basic function o the Prandtl torsional 
stress function ‘’ was used; this satisfies the two-dimensional equation: 


V¥+2 = 0 (2.4,3) 
and vanishes on the closed boundary, so 


Y = —2o. (2.4,4) 
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It was shown to be possible to place very close bounds to the errors in 
approximations to ¥’ itself (1,2) and to a stress function arising in the 
St. Venant theory of flexure (2). 


3. A method for estimating the errors in step-by-step solutions of 
sets of ordinary differential equations 

3.1. Outline of the method 

The method here described appears to have been first proposed by 
L. F. Richardson (3, 4), who gave it as an example of what he called ‘the 
deferred approach to the limit’. A more recent account (5) has been given 
by the present writer and the practical value of the method, which is 
applicable to non-linear equations, has been demonstrated by a number 
of examples. 

Briefly, the method is based on the idea of extrapolation towards the 
limit of the step-by-step solution corresponding to a vanishingly small 
interval. Suppose that the independent variable is ¢, the range of integra- 
tion a to b, and n the number of equal intervals used in the step-by-step 
process. Then, provided that the process is completely regular and the 
boundary values of the unknowns are all given for a (say), we may assume 
that the error in the value of the dependent variable x, can be expanded 
in the series 

e,(t) = n-*(e, +e, n+, n+ etc.), (3.1,1) 


provided also that 7 is not too small. The number k& is characteristic of 
the particular process used and is called its index by the present writer. 
Now when the interval is sufficiently small we may as a first approxima- 
tion retain only the first or dominant term in the expansion and write 


«{t) = en. (3.1,2) 


Since & will be known it becomes possible to calculate e¢, when the approxi- 
mate solution has been obtained for two values of n. Thus we shall have 


a(t) = 2,,(t)+e9nz* = z,,(t)+egng*, 


where 2,,(t), x,,(t) are the approximations to z,(t) with n, and n, intervals 
respectively. The last equations yield 


Cn == 


Xp4(t)—2,,(t) (3.1,3) 
k 


ny *—nz 


and a first approximation to the error is obtained. With three values of 
we can similarly calculate e, and e, and obtain a next approximation, as 
has been shown by examples, and the process could be carried still further 
if desired. 
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It may be noted that the index varies from 1 in the original and relatively 
crude process given by Euler in the eighteenth century to 4 in the process 
of Runge and Kutta (6). The various methods are briefly reviewed in the 
writer’s paper (5). 

When the method just described is applied an estimate should first be 
made of the number of steps needed to provide the desired accuracy. 
Then a first calculation should be made with, say, half this number of 
steps. A comparison of the results obtained with the two numbers of 
steps allows the error to be assessed and partially corrected. If need be, 
a still larger number of intervals must finally be used. 

Evidently the errors introduced by rounding-off the results of a com- 
putation will reduce the accuracy of the results yielded by the foregoing 
process. It will therefore be desirable to minimize such errors by carrying 
the computations some decimal places farther than might appear strictly 
necessary. The greatest error which can be introduced by rounding-off 
can easily be ascertained in any given case. 


3.2. Extension of the method to partial differential equations 

The method just described can sometimes be applied to partial differen- 
tial problems. It is essential that a perfectly regular process be used and 
that an index should exist. When this is so, the error can be assessed as 
before when the results for two similar lattices are known. It may even 
be possible to dispense with this condition of similarity. 


4, Estimation of the errors in approximate solutions of linear 
problems when an approximation to the Green’s function is 
known 
We shall suppose that we have an approximation y%, to the solution of 
2.3.2) which satisfies the linear and homogeneous boundary conditions 

exactly and gives a residual 


e = f—A(¥,+8). (4,1) 


Then, if there is an approximation G,, to the Green’s function appropriate 
to the homogeneous boundary conditions, we may derive an approxima- 
tion to the error of deficiency in %, at any point by forming the corre- 
sponding integral of the product «@, throughout the region. 
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